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LINE-LOADINGS ON FINITE LENGTH CYLINDRICAL 
SHELLS—SOLUTION BY THE FINITE FOURIER 
TRANSFORM+ 


By WILLIAM A. NASH (University of Florida, College of Engineering, 
Gainesville) and 
T. F. BRIDGLAND, Jr. (University of Florida) 
[Received 9 June 1960] 


SUMMARY 


The linear equations due to Fliigge for the thin elastic circular cylindrical shell 
are solved by the use of the finite Fourier transform for two general cases of radial 
line-loadings on a finite length shell. The loadings considered are: (a) an arbitrary 
line-load acting along a generator, and (b) an arbitrary line-load applied around the 
circumference. The solution presents the orthogonal displacement components of 
the middle surface in the form of rapidly converging series, 


1. Introduction 


Tue problem of elastic deformations of a thin circular cylindrical shell 
subject to a radial loading applied along a generator of the shell has been 
examined recently by several investigators. In 1954 Hoff, Kempner, and 
Pohle (1) employed Donnell’s equations (2) to investigate the deforma- 
tions of a simply supported, cylindrical shell of finite length, subject to 
a sinusoidally distributed radial force that is: (a2) symmetric with respect 
to the midpoint of the length of the shell, or (b) antisymmetrically distri- 
buted with respect to this same point. The results of these studies made 
it possible to consider the case of circumferential moments distributed 
along a generator of the cylinder. The problem of a centrally located, 
uniformly distributed radial line load applied along a generator was later 
considered by Kempner (3) on the basis of Fliigge’s equations for circular 
cylindrical shells (4) and Donnell’s equations. The results of solutions of 
the two sets of equations were shown to be in very good agreement for 
the particular geometry investigated. A treatment of this same problem, 
incorporating the effects of transverse shear deformations, was later carried 
out by Cooper (5). The problem of arbitrary radial loads applied along 
a generator of a thin cylindrical shell has been attacked by Sjostrom (6) 
who employed integration by residues to solve the Fliigge equations. 

The problem of deformations of a thin circular cylindrical shell subject 


+ The results were obtained in the course of research sponsored by the Office of Ordnance 
Research, U.S. Army, under Contract DA-01-009-ORD-671 with the University of Florida. 
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to radial line-loads acting circumferentially along a normal section of the 
shell has been treated for the special cases of (a) uniform line-load acting 
on a shell of infinite length, and (b) equal and opposite concentrated radial 
forces. The former study is due to Klosner and Pohle (7) and the latter 
to Yuan and Ting (8). The case of a uniformly distributed segmental 
radial load has been considered by Sheng (9). 

The purpose of the present investigation is to offer a solution to the 
problem of a thin shell of finite length subject to an arbitrary line-load 
applied in a radial direction around the circumference of the shell. The 
circumference corresponds to a section of the shell normal to the geometric 
axis. Further, the present investigation offers a unified approach to the 
problem of the line-loaded cylinder since the shell subject to arbitrary 
radial line-loads along a generator is shown to be amenable to the same 
type of analysis. 


2. Basic equations 


Let us consider a thin circular cylindrical shell of radius a and thickness 
h. Young’s modulus of the material is denoted by E and Poisson’s ratio 
by v. We select a coordinate system in the middle surface such that the x* 
axis coincides with a generator; a® denotes the circumferential coordinate, 
where ® represents the central angle, and the z-direction is radial (taken 
to be positive inward). The corresponding displacement components of 
a point in the middle surface are denoted by u*, v*, and w* respectively. 
It is convenient to work with dimensionless displacement components u, 
v, and w, defined by the relations u = u*/a, v = v*/a, and w = w*/a, as 
well as with the dimensionless coordinates x = x*/a and ®. 

Fliigge’s differential equations of equilibrium for circular cylindrical 
shells (4) may be written in the dimensionless form 


2 


where dots indicate differentiation with respect to ® and primes denote 
differentiation with respect to x. Here q is the applied surface load acting 
normal to the middle surface and taken to be positive in the inward radial 
direction. Also, k®? = h?/12a? and D = Eh?/12(1—v*). 
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The finite Fourier transform of a function f(t) is defined by the relation 
(10) 


Ff} = = srk). (4) 
The inverse transform is defined ia 
= =f. (5) 
Also, from (10) we have the relation 
= ikf *(k) (6) 
as well as analogous relations involving higher derivatives. 


3. Radiai forces acting along a generator 
Let us represent the normal load in the form 


= a,sin("™)) (7 


n=1 


and the displacement components by the relations 


w(x,) = F,,(@)sin( 7"), (8) 
v(x, 0) = > (9) 


n=1 
Henceforth we write A = nza/L and, for the case of a radial line-load 
acting along a generator ®,, we set 


f@)=1 for (11) 
For such a loading we have from (4) 
1 +e 
ge — £o-ire,(Sin re 


The radial load, corresponding to (11) and (12), is P = 2ae, hence 


$40) = (14) 


| 

? 

b 
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and from (5) f(®) = (15) 


re 
as « > 0, (15) becomes 


f@) = >, (16) 


We denote the finite Fourier transforms of the functions F,,,, appearing 
in (8), (9), and (10) by 
nn (®)} = Fin(8) 
and further set P, = Pa, and also take is = p. The transforms of (1), 
(2), and (3) thus become 


pF) = 0, 


Fi, (s) = 0. 
The solution of this set of linear equations will be denoted by 
Finals) = An: (20) 


m=1 


Thus, for a loading 


(21) 


where (22) 


g(x) being the load per unit length along the generator, we have displace- 
ment components uw, v, w given by equations (8), (9), and (10) where 


Fan(®) 


Finn(8) = Amn()/A,(8). 


Lin 

(23) 

and (24) 
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The terms in (24) are found from (17), (18), and (19) to be 


v*)a, 


—isDo 
l—vp 2r7Eh 


x | 24-201 (T 7) (143k2)n!, 


l—v 27 Eh 


l—y 2a Eh L 


2A, 


Thus we have Fs = (29) 


where H*.,,(s) is obtained by taking quotients of the appropriate quantities 

under the vinculum designations in equations (25) to (28). For m = 1 

and 3, H%,,(s) is an even, real valued function of s whereas for m = 2, 
FE}, (8) is an odd, purely imaginary function of s. Thus we may write 


F,,(0) = 


F,,(®) = TER? (21 (31) 


4 
(26) 
? 
| 4 
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from which the characteristics of the deformations due to the line-load 
along the generator ® = ®, are apparent. 


4. Radial forces acting along a circumferential band 


To treat line-loads acting radially along a circumferential band of the 
shell, i.e. along a section normal to the geometric axes, it is only necessary 
to replace in equations (25) to (27) a, e-*®/2a by 


2a 
where f*(s) = do, (32) 
0 


where f(®) is the applied load per unit length along the section x = z, 
of the cylinder. Thus, for band loading the expression analogous to (29) is 


so that from (23) we have 
2 


If, for example, f *(s) is an even, real valued function of s, then f*(s) Z*,,,(s) 
is even and real for m = 1 and 3; odd and pure imaginary for m = 2. 
In this case, since the displacements u, v, w must be real, we have 


2(1—v®)a?/ . 
F,,(0) = 7 *(0) 0) +2 


The displacement components are now found by substituting (35) and 
(36) in equations (8) to (10). Analogous expressions may be derived for 
the case when f*(s) is purely imaginary. In the speciai «.cisymmetric case 
of a uniform load around the normal section it is convenient to write 
f*(8) P*%,,, 


where P* represents the force per unit length and 4,, is the Kronecker 
delta. 


5. Numerical example 


The problem of a simply supported cylinder free from applied loads 
except for a centrally located, uniformly distributed, inward radial line 
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load over a segment of the generator ®, was considered. The parameters 
were taken to be 

@ = 35 in. h = 0-625 in. L = 140 in. 

E = 30x 10 Ib/in® v= 03 25 = 7 in. 
where 25 is the distance over which the line load is distributed. If P 
represents the resultant load, then from equations (30), (31), (8), (9), and 
(10) the radial displacements at the midpoint of the length of the shell 


and under the load, diametrically opposite the load, and midway between 
these two points, are given by 


ERE 2, Ga 40 |x 


 (—1) Bhan (38) 
respectively, where = 


Bt,(s) = 88+ (2-4674n2—2)s® + (2-2830n*—4-5647n2-+ 1)s4+ 
+-(0-9388n4 — 2-2830n?-+- 2-0973)n2s? 

+ (0-1448n4*—0-1408n2-+ 13030-3315)n. (41) 
The values of (w/P) x 10°, given by equations (37), (38), and (39), are 
found to be 0-158, —0-008, and 0-010 respectively. These values were 
established through the use of the first ten integral values of n and s. 
The use of additional terms did not alter these values in the third decimal 
place. The first of these values, corresponding to the radial deflexions 
under the midpoint of the line load, is to be compared to the value of 
0-159 found by Kempner (3) for the same geometric and elastic parameters. 
The technique employed in (3) consists in considering the displacement 
to be composed of the displacements of two fictitious sheets winding 
around the cylinder an infinite number of times, one clockwise and the 
other anticlockwise. Usually, the first term, corresponding to the first 
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turn of one of these sheets around the shell, suffices for a satisfactory 
approximation to the deflexion. 
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THE SOUTHWELL METHOD FOR PREDICTING 
CRITICAL LOADS OF ELASTIC STRUCTURES 


By 8. T. ARIARATNAM 
(Nunavil West, Chavakachchery, Ceylon) 


[Received 27 August 1959. Revise received 30 April 1960) 


SUMMARY 


The experimental determination of the critical loads of isolated pin-ended struts 
by the Southwell plot is well known. In this paper it is shown that the Southwell 
method may be extended to the prediction of the elastic critical loads of plane 
frameworks. 


1. Introduction 


A critica. load for an elastic structure is one at which a straight form and 
an adjacent distorted form of equilibrium are both possible (small deflexion 
theory). In general, a structure has an infinite number of critical loads but 
it is the smallest of these that is of practical importance. Theoretically, 
therefore, if the loading on an elastic structure is increased slowly, in a 
constant proportion, the structure will remain straight and undistorted 
till the smallest critical load is attained, at which it is capable of assuming 
a deflected configuration of (neutral) equilibrium. In practice, however, on 
account of unavoidable imperfections of manufacture, the structure will 
begin to deflect from the first application of the load, the deflexions and 
stresses becoming large as the smallest critical load is approached. The 
theoretical critical load, based on elastic considerations, is therefore never 
attained as the stresses in the structure would have exceeded the limit of 
proportionality before the attainment of this load. 

Southwell (1) showed that the load-deflexion curve of an initially im- 
perfect pin-ended strut is a hyperbola in the neighbourhood of the smallest 
critical load P., the asymptote of the hyperbola being the horizontal line 
P = P.. By asuitable transformation of variables, this hyperbolic portion 
of the load-deflexion curve may be converted into a straight line whose 
gradient is the reciprocal of the smallest critical load of the strut. In this 
way, it is possible to predict the value of the critical load, although it 
may not be experimentally attainable. 

The Southwell method has recently been employed by a number of 
investigators for the prediction of the smallest elastic critical load of plane 
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frameworks loaded in their plane (2, 3,4). This paper provides a theoretical 
justification of the validity of this procedure not only for buckling of 
structures within their plane but also for the case of lateral buckling out 
of their plane. 


2. Flexural buckling of isolated struts 
Southwell (1) presented his analysis for the special case of a pin-ended 
strut of constant flexural rigidity. In this section it is shown that South- 
well’s analysis may be extended to the more general case of a strut with 
varying flexural rigidity and different combinations of end conditions. 
P x P 


Consider a strut of varying flexural rigidity EJ,, loaded concentrically 
by an axial compressive force P. With axes chosen as in Fig. 1, the equa- 
tion of flexure of an initially straight strut is 

(El, = 0, (2.1) 


where y is the transverse deflexion and primes denote differentiation with 
respect tox. Writing EJ, = where is a differentiable function 
of x, and denoting the ratio P/ EI by A, this equation may be written 


[Viw)y"]"+Ay” = 0. (2.2) 


The general solution of this differential equation is 


Y = A) 04, (2.3) 


where a,..., a, are real constants determined from the conditions of end 
restraint, and ¢,(7,A), ¢,(z,A) are transcendental functions of x and 
the parameter A. To evaluate a,,..., a, we consider the following end 
conditions: 

pinnedend y=0O and y" = 0 

fixed end y=0 and y =0 > (2.4) 

free end y"=0 and = 0 


For any given case, if we substitute the solution (2.3) in the relevant four 
end conditions, we obtain a set of four linear, homogeneous equations in 
the constants «. Non-trivial solutions of this set of equations will exist, 
indicating existence of deflected configurations of equilibrium if and only 
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if the determinant A(A) of the coefficient matrix vanishes. In general, 
the equation A(A) = 0, being a transcendental equation in A, gives an 
infinite number of roots A, (r = 1, 2,...) defining an infinite number of 
critical loads P. = EIA,. Associated with each root A, is a set of ratios 
= by (x, A,) + + Say] (2.5) 
indeterminate as to magnitude but defined as to shape. The critical loads 
P. are thus defined in terms of the characteristic values A, of an eigenvalue 
problem whose characteristic functions are the buckling modes y,. 
2.1. Orthogonality of the buckling modes 
The characteristic function y,(x) satisfies the differential equation 
= 0. (2.6) 
Multiplying this equation by the characteristic function y,(z) and integrat- 


ing over the length L of the strut, we obtain 
L 


de = 0. (2.7) 
0 
Integrating by parts, this becomes 
L L 
0 0 


The right-hand side vanishes for any combination of the linear homo- 
geneous boundary conditions (2.4), and therefore 


L 
dx =A, f de. (2.8) 
0 


Interchanging y, and gh we obtain 


dx = af dz. (2.9) 
0 


Subtraction of one equation from the other gives 


(A,— A,) yy, dx = 0. 
Ifr AsandaA,#A,, then 


de = 0. (2.10) 
é 


We deduce immediately from equation (2.8) that 


L 
de = 0. (2.11) 
0 


4 
m 
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When r = s, we obtain from equation (2.8) the normalizing relation 


L L 
| de dz. (2.12) 
0 


For convenience we can, by suitable choice of the constant a,,, make each 
side of this equation equal to A,/ ZI. 

Equations (2.10) and (2.11) are the orthogonality relations satisfied by 
the buckling modes y,. 

It can be shown (5) that the set of all characteristics functions y, of 
the differential equation (2.2) for a given set of boundary conditions of 
the form (2.4) comprise a complete system. We may therefore expand any 
arbitrary deflexion function y(x) satisfying the same geometric boundary 
conditions as the functions y, as an absolutely and uniformly convergent 
series (6) of the form a 
= 2 


where a, are real constants given by 


L 
a, = EI { 
0 
2.2. Magnification of imperfections 
Consider a strut with an initial imperfection y,(z) which, from the 
previous section, may be expressed as an infinite series in terms of its 
characteristic deflexion modes y,(x) in the form 


= a,y,(2. (2.13) 


When the strut is loaded by an axial compressive force P, let its 
deflexion from the straight position be y(x) given by 


y(z) = blz), (2.14) 


where the b, are real constants to be determined. 
To obtain 6, we use the theorem of minimum potential energy (cf. 7). 
The total potential energy of the strut under an axial force P is given by 
L L 
U=4 Eldy’—y)* de—4 de. 
0 0 
Using the relations (2.10), (2.11), and (2.12) after substitution from 
(2.13) and (2.14), this can be written as 


1 +-constant. 
1 
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The change in the potential energy when the strut is given a virtual 
displacement 5y = > 8, y, from its equilibrium position is 


sU = 2 (8b,)*(A,—A). 
By the theorem of stationary potential energy, the condition for equi- 


A,(b,—a,)—b, = 0, 
a a 

ich gi =—+ 2.15 
which gives b, I-A, ~ 1- PP for all r ( ) 
The equilibrium is stable if and only if this stationary value for U is a 


minimum, i.e. if - ; 
2 (5b,)*(A,—A) > 0, 
or A<A, for all r. (2.16) i 
Substituting from (2.15) into (2.14), the transverse deflexion of the strut 
under an axial thrust P is given by 


a, 


ue) = > (2.17) | 

2.3. The Southwell plot 
The deflexion of the strut measured from the initially curved position is 


a, 


If the experimental deflexion A is measured at the point x = 2, (usually 
the centre), 


a, 
(2.19) 


In practice (where y, arises from imperfections of manufacture), we 
may expect the first term a, y,(x,) to be of the same order of magnitude 
AS @,Y(Xo), ete. Hence when P is comparable with the smallest critical 
load P,, and P, > P,, we may write 


(2.20) 
1 


which is the equation of a rectangular hyperbola in the variables P and 
A. Denoting A/P by v, this equation may be transformed to 


_ A, ‘ 
bed == (2.21) 
This is a linear equation connecting v and A, so that when observations 
of A and P are plotted on a graph relating v and A, they should lie on 


ig 
3 
q 
AP 
or 
ee 
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a straight line which cuts the axis v = 0 at the point A = —a, y,(z,) and 
has gradient equal to the reciprocal of P,. This line is called the ‘South- 
well line’. 


3. Flexural buckling of plane frameworks 

We consider plane, rigid-jointed, frameworks (Fig. 2) made up of 
straight members of uniform cross-section, loaded in their plane by ex- 
ternal loads applied at the joints so as to 
produce axial forces only and no bending 
moments in the members. We assume that 
the directions of the external loads are not 
altered by any deflexion of the framework, 
and that the individual loads bear a con- 
stant ratio to one another as the general 
intensity is increased, i.e. the loading is 
proportional. Let A denote the load para- 
meter. Further, we suppose that the frame- 
work is sufficiently restrained against 
lateral buckling so that all the deflexion of the framework is confined to 
its plane. 

With these assumptions and on the basis of a linearized small deflexion 
theory it can be shown [8 (a)] that at certain critical values A, (r = 1, 2...., 
A... < A, say) of the load parameter A, given by the solutions of a trans- 
cendental equation A(A) = 0, the framework can assume a neighbouring 
deflected configuration of neutral equilibrium. We call A, the rth critical 
load parameter of the framework for the given loading pattern, and its 
associated deflected shape »,, the rth buckling mode. The function », is 
considered to denote collectively the transverse deflexion functions of all 
the members of the framework. The buckling modes », are obtained from 
the non-trivial solutions of a set of linear homogeneous equations and are 
therefore indeterminate in magnitude but defined as to shape. 


3.1. Orthogonality of the buckling modes 
The buckling mode », corresponding to the critical load parameter A, 
satisfies the differential equation 
(EIn;)" +A, Py, = 0 (3.1) 
for every member. Here FJ is the flexural rigidity of a typical member 
of length L, AP its axial compressive force at load parameter A, and primes 
denote differentiation with respect to z, where z is measured from one 
end along the member (Fig. 3). 
Multiplying equation (3.1) by »,, where 7, is the buckling mode 


MP, AP, 


J 
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associated with the different critical load parameter A,, and integrating 
over the length of the member, we obtain 


L 
+A, dx = 0. 
0 


Twofold integration by parts of the first term and single integration by 
parts of the second term give 
L 


0 
(M, 9.)2-0+ (M, (S, A,)z-0+ (S, A,)2=L 


(3.2) 
A B 
Ae 


Fic. 3 


where the quantities M, S, @, A are shown in Fig. 3. Summation of 
equation (3.2) over all the members of the framework gives 


L L 
0 


0 


The right-hand side vanishes on using the boundary conditions (2.4) in 
conjunction with the principle of virtual work.t We have, therefore, 


2,3 / Pn, de. (3.4) 
Interchanging r and s, we get 
L 
= Pr), x, dr. (3.5) 
Subtraction of one equation from the other leads to 
(A,—A,) Pn, 7, dx = 0. 


+ It is assumed that the framework is so supported that the reactions at the external 
supports do no work when the structure deforms. 
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a straight line which cuts the axis v = 0 at the point A = —a, y,(z) and 
has gradient equal to the reciprocal of P,. This line is called the ‘South- 
well line’. 


3. Flexural buckling of plane frameworks 

We consider plane, rigid-jointed, frameworks (Fig. 2) made up of 
straight members of uniform cross-section, loaded in their plane by ex- 
ternal loads applied at the joints so as to 
produce axial forces only and no bending 
moments in the members. We assume that 
the directions of the external loads are not 
altered by any deflexion of the framework, 
and that the individual loads bear a con- 
stant ratio to one another as the general 
intensity is increased, i.e. the loading is 
proportional. Let A denote the load para- 
meter. Further, we suppose that the frame- 
work is sufficiently restrained against 
lateral buckling so that all the deflexion of the framework is confined to 
its plane. 

With these assumptions and on the basis of a linearized small deflexion 
theory it can be shown [8 (a)] that at certain critical values A, (r = 1, 2...., 
A,-1 < A, say) of the load parameter A, given by the solutions of a trans- 
cendental equation A(A) = 0, the framework can assume a neighbouring 
deflected configuration of neutral equilibrium. We call A, the rth critical 
load parameter of the framework for the given loading pattern, and its 
associated deflected shape »,, the rth buckling mode. The function 7, is 
considered to denote collectively the transverse deflexion functions of all 
the members of the framework. The buckling modes », are obtained from 
the non-trivial solutions of a set of linear homogeneous equations and are 
therefore indeterminate in magnitude but defined as to shape. 


3.1. Orthogonality of the buckling modes 
The buckling mode », corresponding to the critical load parameter A, 
satisfies the differential equation 


(EIn;)"+A, Pn, = 0 (3.1) 


for every member. Here £J is the flexural rigidity of a typical member 
of length L, AP its axial compressive force at load parameter A, and primes 
denote differentiation with respect to x, where x is measured from one 
end along the member (Fig. 3). 

Multiplying equation (3.1) by »,, where », is the buckling mode 


| 
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associated with the different critical load parameter A,, and integrating 
over the length of the member, we obtain 


L 
+A, Pot dx = 0. 
0 


Twofold integration by parts of the first term and single integration by 
parts of the second term give 
L 


Pat ni) de = (EI +4, 
0 


(3.2) 


Fic. 3 


where the quantities M, S, 6, A are shown in Fig. 3. Summation of 
— (3.2) over all the members of the framework gives 


The right-hand side vanishes on using the boundary conditions (2.4) in 
conjunction with the — of virtual work.t We have, therefore, 


Interchanging r and s, we get 


Subtraction of one equation from the other leads to 


| 


+ It is assumed that the framework is so supported that the reactions at the external 
supports do no work when the structure deforms. 
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If A, # A,, we have the relation 


Pn, 7, dx 


and, from equation (3.4), the further relation 


L 
dx = 0. 
0 
When r = 8, A, = A, and we obtain the normalizing condition 


L L 
0 0 


in which, by a suitable choice of the magnitude of n,, each side may be 
made equal to A,. 

It may be shown (9) that the buckling modes », corresponding to all 
the critical load parameters A, (r = 1, 2,...) form a complete set of ortho- 
gonal functions, so that we may expand any arbitrary deflected shape y 
of the framework (satisfying the same geometric boundary conditions as 
the functions »,) as an absolutely and uniformly convergent series of the 
form 

y = 2% Np» 


where a, are real constants given by 
L 

a, = > Py'n, dx, 


> denoting summation over all the members of the framework. 


3.2. Magnification of imperfections 
Consider a rigid-jointed framework whose members have an initial 
deflexion from the straight position (due to imperfections of manufacture) 


represented by the deflexion function y). We may express the function 
Yo a8 a series of the form 


Yo = Vrs (3.10) 


where », is the buckling mode corresponding to the rth critical load para- 
meter A,, and the real constants a, are obtainable from relation (3.9). 

When the framework is loaded (as in Fig. 2), let the deflexions from 
the initially straight positions of the members, at a value A of the load 
parameter, be given by the function 


¥= 26% (3.11) 
where , are real constants to be determined. 


=0 (3.6) 
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The total potential energy of the framework in its deflected configura- 
tion is 


L 
U=> dx-+constant, 


if we ignore the effect of direct axial strain. Substituting the expansions 


(3.10) and (3.11), and using the relations (3.6), (3.7), and (3.8), this 
becomes 


U = 


Consider a virtual displacement of the framework from its equilibrium 
position, represented by ‘ia 
dy = 2%, 


The change in the potential energy due to this virtual displacement is 
BU = (6b,)*A,—d). 

By the theorem of stationary potential energy, the condition for equi- 

librium is A,(b,—a,)—b, = 0 

which gives b, = a,/(1—A/A,) for all r, (3.12) 

and, by the theorem of minimum potential energy, the equilibrium is 


stable if and only if A < A, for all r. 
Substituting from (3.12) in (3.10) the equilibrium deflexion y is there- 


fore given by 
y= —— Nr 
he 1—A/A, 


It follows from this result and by the discussion in section 2.3, that, 
provided a, # 0 and A, > A,, the Southwell plot may be employed for 
the prediction of the first critical load parameter A, from the experimental 
observations of external load and transverse deflexion at any suitable 
point of a general member of the framework. 


4. Torsional-fiexural buckling of plane frameworks 


In frameworks whose members are deep in the plane of the framework 
and shallow in the perpendicular direction, buckling of the framework 
out of its plane may precede buckling within its plane. In the special case 
where the buckling modes do not involve any lateral displacement of the 
joints, the stability equations for determining the critical loads may be 
set up in the manner indicated by Masur and Cukurs (10). If the joints 
are not restrained against lateral movement, the smallest critical load is 
the one associated with joint translation and the stability equations may 
be then set up as indicated below. Here again the external loading is 
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considered to consist of point loads acting at the joints of the framework 
and in its plane so that no bending or twisting moments are induced in 
the members before buckling occurs. Further, it is assumed that the 
framework is supported in such a manner that the external reaction forces 
will do no work when the framework deforms. 

We shall establish the validity of the Southwell plot in the case when 
the members of the framework are of bisymmetrical cross-section and 
possess large warping rigidity, e.g. rolled steel joists. When the members 
are of monosymmetrical cross-section, such as double angle or T-sections, 
the proof is similar. 


4.1. Equations for a single member 
Let a typical member of length L between joints i and j buckle laterally 
as shown in Fig. 4. The displacements (u,v) of the shear centre (which in 
this case is coincident with the centroid) and the rota- 
tion 8 of the cross-section are given by the differential 

equations [cf. 8(5)], 
El,u*+APu" = 0, (4.1 a) 
EI, v'¥+APv" = 0, (4.1 b) 
ETB +(APp?—GK)p” = 0. (4.1 ¢) 
buckled cross- In these equations AP is the axial compressive force in 
—— the member along its centroidal axis when the load 
parameter of the external loading is A. EI,, EI, represent the relevant 
flexural rigidities, and ET, GK are the warping and torsional rigidities of 
the cross-section respectively; p is the radius of gyration of the cross- 


x-y plane (plane of framework) plane 
Fic. 5 


section about the shear centre, and primes denote differentiation with 
respect to the coordinate x which is measured from one end along the 
member (Fig. 5). 

The buckling in the plane of the framework (the zy-plane) is uncoupled 
as it is clearly possible to have pure bending in the zy-plane unaccom- 


yw i x j 
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panied by any twist or lateral deflexions in the members. We may there- 
fore ignore equation (4.1a). The solutions of equations (4.1 b) and (4.1 ¢) 
may be expressed in slope-deflexion form. 

Denoting by @,, 6,, z the angle of twist, bending rotation and lateral 
displacement of the ends, the twisting moment M,, the bending moment 
M,, and the shearing force Q, at the two ends i and j are given by 


M,, = 88,,—88,, 

and M,; = —S0,,+58,, 


6EI 6EI 
Q4 = $1) LB 
where 4,,..., 6, are the stability functions due to Livesley (11) and the 
torsional stiffness S (see Appendix) is given by 


p(GK —APp?) GK—APp* 


Using matrix notation, we may rewrite equations (4.2) and (4.3) in the 


S= 


form F, == Vu D,+ Yi; D; (4.4) 
F; = Y;,D,+ ¥,,;D, 
M, 6, 
Q. 
and Yu = 0 0 
4EI 6EI 
4 
6EI 
Y,=f[-s 0 6. 
6EI 
0 L $s 
6EI 12z1 


f 
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0 


4.2. Stability equations for a framework 

Equations (4.4) are in terms of axes (x,y,z) chosen with respect to a 
member. To build up the stability equations for a frame, we must rewrite 
them in terms of a fixed set of axes (x’, y’,z) (called frame axes) common 
to all members of the framework. 

Let the inclination of the member ij to the x’ axis be « (see Fig. 5). 
If F’, D’ with relevant subscripts denote the force and displacement 
systems at the ends of the member in the primed system of coordinates 
(viz.: 2’, y’,z), then, 


F,= TF; F,=TF, and D,= TD; D,=TD; 


where T is the orthogonal transformation matrix 


cosa sina 0 
—sina cosa 
0 0 1 
Equations (4.4) in frame coordinates then become 
Fi = 
Fj = ¥j,Di+¥j,D; 
where Y}, = T*Y,, T (r, s = i,j), and T* 
denotes the transpose of T. 

Consider a typical joint i at which three 
members a, 6, ¢ are rigidly connected 
together, Fig. 6. If A; is the displacement 
vector of joint i (in frame coordinates), 
the condition of compatibility of displace- 
ments gives 

(Di)q = (Di), = (Di). = (4.6) 
Since the external force system consists of concentrated loads acting at 


(4.5) 
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the joints in the plane of the framework, the condition for equilibrium 
of the force system at i gives 


(Fiat+(Fint (File = 9. (4.7) 
Substituting from equations (4.5) in (4.7) we get 
Vis Ay + Vu t+ = 0. (4.8) 
If the framework has n joints, we can write n sets of equations similar 
to (4.8) giving a total of 3n scalar equations. We may write these equa- 
tions in the form YA = 0, (4.9) 
where A is the vector A, 


A, 


A, 
and Y is the symmetrical square matrix obtained by combining the stiff- 
ness matrices of the individual members as in equation (4.8). 

Equations (4.9) form a linear homogeneous system possessing non- 
trivial solutions if and only if det Y = 0. det Y is a transcendental function 
of A and this condition gives, in general, an infinite number of roots A, 
(r = 1,2,...; A... <A,, say) which are the critical load parameters of the 
framework at which lateral buckling will occur. Associated with each of 
these values A, is a transverse deflexion function v,(z) and a rotation 
function £,(x), both defined over the lengths of the members. A, is the 
rth critical load parameter of the framework and the deflected configura- 
tion defined by the functions v,(x), B(x) is the rth buckling mode. 


4.3. Orthogonality of the buckling modes 


The functions v,(x), 8,(x) corresponding to the critical load parameter 
A, satisfy the differential equations 


EI, vi’ +A, Pos = 0, 
+(A, = 0, 
for each member. 
Multiplying the first equation by v,(x) and the second by £,(x), where 
(v,,8,) correspond to a different critical load parameter A,, and integrating 


both equations over the length L of the member, we obtain 
L 


(BI, +A, dx 


0 


L 
= [(El, v7 +A, + (El, dx 
= 0, 0 


* 
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4 
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L 
| +, dx = [{ET +(A, — 


L 
+ | (ETB; Pp*—GK)B, 
== 
Adding these two equations and summing over all the members of the 
framework give 
L 
> (Bh GKB, B,—A, Plv,v, + dz 


+ [{(@K—A, Pp*)8,— ET + [ETB 

By the principle of virtual work, the sum of the first three terms on the 

right-hand side is zero. The fourth term is also zero if the end cross- 

sections of the members do not warp, so that fp’ = 0 at the ends (see 
appendix). The above equation therefore simplifies to 


L L 
> | (BL GKB, B,) dex = 2, (4.10) 
0 0 


Interchanging r and ¢ in equation (4.10) leads to 


L L 
0 0 


Subtracting equation (4.11) from equation (4.10), we obtain 
L 
0 


If r #8, A, #A, and we have the orthogonality relation 


L 
+8, = 0, (4.12) 
0 


which on substitution in equation (4.10) or (4.11) gives the additional 
relation 


L 
(EI, ETB, B,) dx = 0. (4.13) 
0 


When r = 3, A, = A, and we obtain from equation (4.10) the normalizing 
relation 


L L 
| (Blo? + + OKB?) de = A, [ (4.14) 
0 0 
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each side of which, by suitable choice of the magnitude of the functions 
v,(x), B(x), may be made equal to A,. 

If v(x) is a given deflexion function of the members of the framework 
and §(2x) is the associated rotation function of the members satisfying the 
compatibility conditions (4.6) at every joint, we may expand these func- 
tions as series in terms of the characteristic functions v,(x), 8,(x) in the 


formt x 
v(x) v,(x ), B(x )= 2 B(x), 


where a, are real constants given by 


L 
a, = Plo’ v,+p°6'B,) de, (4.15) 
0 


taking account of the orthogonality condition (4.12) and the normalizing 
relation (4.14). 


4.4. Magnification of imperfections 
Suppose that the members of the framework are initially not straight 


and have small imperfections given by the functions v(x), By(z). We may 
express these functions in the form 


= a, v(x), Bo(x) = a, Bx), (4.16) 


where the real constants a, are given by (4.15). 

On application of external loads, given by load parameter A, let the 
deflected configuration of equilibrium of the framework be given by 
functions v(x), B(x) measured from the initial straight position, where 


= Bl) = (b, real constants). (4.17) 


The total potential energy of the framework is 


L 
+-constant, 
+ Note that the corresponding oa in the expansions for v(x) and B(x) may be 
taken as equal; for if we take f(z) = 2 Pr Be), then, since the boundary conditions 


(4.6) are linear and homogeneous, the functions v(z) = 0 and A(x) = > (bp —@,) B(x) also 
r=1 


satisfy the compatibility relations (4.6). But if v = 0, by (4.6), 8 = 0 at the ends of the 
members (except in the trivial case of a straight continuous beam) and therefore B is 
identically zero everywhere. Hence 


=0, giving 2 br Bla) = Brie). 
rT = 


> 
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where 5 denotes summation over all the members of the framework. 


m 
Substituting from (4.16) and (4.17) and using the relations (4.12), (4.13), 
and (4.14) this becomes 


U = constant. 
r=i 
Subjecting the framework to a virtual displacement = 5b, v,, 
58 = > 5b,8,, the potential energy changes by an amount 
r=1 


The theorem of minimum potential energy then gives 
b, = a,(1—aA/A,)“1, for all r, (4.18) 
for equilibrium, and A < A,, for all r, for stability. 
Substituting from (4.18) in (4.17), we obtain 


@ 


a, 


oe) = Ble) = 


showing that the Southwell method may be used for the experimental 
prediction of the lateral buckling loads also. 
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APPENDIX 
Solution of the Equation of Torsion 


The torsion equation (4.1 ¢) is 
+ (APp*— GK)p” = 0 

and has the solution 

B = Acoshuxr+ Bsinhur+Cr+D, 
where A, B, C, D are real constants and nu? = (@K—APp*)/ET. The end conditions 
B=B,atz=0, and atz= L. 
If the members are rigidly joined together at their ends, we may assume that the 
end cross-sections do not warp, so that 8’ = 0 at the ends. 

Using these conditions to evaluate the constants A, B, C, D, we obtain 
M,, = (@K—APp*)p'(0)— ETB”"(0) = — ETp*B 

= p(GK —APp*)(B,—B;)(2 tanh pL — pl) 


= = 
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ELASTIC WAVES IN A NATURALLY CURVED RODt 
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SUMMARY 

The propagation is considered of free elastic waves of small amplitude in a naturally 
curved rod where the neutral axis forms a plane curve of constant radius. The 
effects of rotatory inertia and radial shear are included, in the same way as in 
Timoshenko’s theory for straight rods, as well as the effect of the extension of the 
neutral axis. 

When the curvature is slight there are considerable simplifications and a Timo- 
shenko-type equation is obtained which governs the flexural motion. A simple 
relationship is then found with the phase velocity of flexural waves in the straight 
rod. This means that data already available for straight rods are now readily 
applicable to slightly curved rods of the same material and cross-section. 


1. Introduction 

A TrmosHENKO-type theory is here used to study the propagation of free 
elastic waves of small amplitude in a naturally curved rod where the 
neutral axis forms a plane curve of constant radius. The uniform cross- 
section is symmetrical about this plane and there is no motion perpen- 


dicular to it. Such rods are frequently used in engineering structures and 
a knowledge of their behaviour is also of value when considering the more 
complex problems of wave propagation in circular cylinders. The analysis 
was stimulated by interest in problems arising from the behaviour of 
structures just subsequent to static fracture and where subsidiary damage 
was encountered as a result of dispersion of the ensuing stress waves. 
In view of the possible more general interest of the present analysis, 
however, it is presented independently of discussion of these particular 
problems. 

In the elementary theory, e.g. Love (1), the effects of rotatory inertia 
and radial shear are neglected. In addition, it is assumed that the central 
line (i.e. the line joining the centroids) remains unextended during the 
flexural motion, while the flexural characteristics are ignored when con- 
sidering the extensional behaviour. These assumptions are often satis- 
factory for low-frequency modes of vibration in rods which are only 
slightly curved, but they are not necessarily valid for the high-frequency 
wave motions such as are initiated by a sudden disturbance. 

Waltking (2) sought to advance the elementary theory by considering 
the effect of the extension of the central line on the flexural motion, while 
om i Crown Copyright reserved, reproduced by permission of the Controller, H.M. Stationery 
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Philipson (3) took into account also the rotatory inertia. However, ex- 
perience in the theory of straight rods makes it clear that the correction 
for radial shear cannot be neglected. 

In the present paper, a unified theory is given which includes the effects 
of rotatory inertia and radial shear in the same way as Timoshenko’s 
(4, 5) theory for straight rods and includes also the effect of the extension 
of the neutral axis. The concepts of standard shell theory, e.g. Herrmann 
and Mirsky (6), are followed and it is shown that when the curvature is 
slight it is possible to make considerable simplifications to the equations 
and that the extension of the neutral axis has no effect upon the flexural 
motion. Thus a Timoshenko-type equation is obtained to govern the 
flexural motion. It is then found that there is a simple relationship with 
the phase velocity of flexural waves in the straight rod of the same 
material and cross-section. This relationship is just the same as that 
predicted by the elementary theories for straight and slightly curved rods 
and it is important because it means that the data already available for 
straight rods are now readily applicable to the slightly curved rod. 

In the case of pronounced curvature, the accuracy requires investiga- 
tion by making comparison with experiment or exact theory. However, 
the only data available are given in an exact analysis by Gazis (7) for 
the plane strain vibrations of a thick-walled circular cylinder but, un- 
fortunately for the present purpose, Gazis’s numerical results are limited 
to a low-frequency mode and so it is possible to make only a very limited 
assessment. Notwithstanding, it is found that for this low frequency the 
present approximate results are less than 5 per cent in error even for a 
cylinder with a wall-thickness equal to the internal radius. 

With regard to the extensional waves it is well known in the case of 
straight rods, e.g. Abramson et al. (8), that the restrictive assumptions 
such as are imposed here, and in the elementary theory, upon the deforma- 
tion of the cross-section do not permit a satisfactory investigation except 
when the wavelength is long compared with the dimensions of the cross- 
section. Within this limitation, it is found that the velocity of extensional 
waves along the neutral axis is predicted sufficiently accurately by 
elementary theory. 


2. Notation 


internal radius of circular cylinder under plane strain. 
cross-sectional area of rod. 

external radius of circular cylinder under plane strain. 
velocity of extensional stress waves in a straight rod 
according to elementary theory, c = (Z/p)*. 


3 q 
= 
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phase velocity of wave propagation along the neutral axis. 


(G/p)t. 
E Young’s modulus of elasticity. of 
Ep increased modulus, Ep, = E/(1—v?). 
G shear modulus. 
h wall-thickness of cylinder under plane strain, h = b—a. 

k geometrical constant defined by equation (11). 

ky, ky geometrical constants defined by equations (13). 

K, K*, K** geometrical constants defined by equations (18). 

L wavelength measured along the neutral axis. 

M ' yesultant moment referred to the neutral axis. 

N resultant circumferential force. 

Pp natural frequency. 

Q resultant radial shearing force. 

R radius of neutral axis. 


non-dimensional radial coordinate, r = r’/ R. 

circumferential coordinate measured around the neutral 

axic. 

t time. 

v, Ww non-dimensional displacements at the neutral axis, see 
equation (8). 

non-dimensional displacements, = ® = @'/ R. 


~ 


components of displacement respectively in the s, z divec- 

tions. 
z non-dimensional coordinate, z = z’/ R. bd 
2" coordinate measured along the inward drawn normal. 
€s> Yre circumferential and shearing strains respectively. 
¢ change in slope of the normal to the neutral axis. 
shear constant. 
A wave number. 
v Poisson’s ratio. 
p density. 
Ta Sys Tre plane polar stress components. 


3. General considerations 

The naturally curved rod is shown in Fig. 1. It is assumed that the 
neutral axis is a plane curve of radius R. The circumferential coordinate 
measured around the neutral axis is s’, while the coordinate z’ is measured 
from the neutral axis along the inward pointing normal and is related 
to the general radial coordinate r’ by the expression 


r’ = 


cp 
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It is convenient to introduce a non-dimensional system of coordinates 
defined by the equations 

== r=r'/R, s=a'/R, (1) 
so that r= l—z. (2) 
The magnitude of the radius R of the neutral axis is determined, as in the 
usual static theory (9), by the equation 
0, (3) 
‘4 


so that as the curvature increases the neutral axis moves farther away 
from the central line and towards the centre of curvature. In equation 


6 


N.8. r=r/R 
s=s/R 
z=2J/R 


Area A 


Neutra! oxis 


Fie. 1. The naturally curved rod. 


(3), the integral is taken over the whole area A of the uniform cross- 
section. 

The non-dimensional displacements in the s- and z-directions are de- 
noted respectively by 6 and #, where 

= 0'/R, = @'/R. (4) 

The equations of motion of the elemental slice shown in Fig. 2 are obtained 
by taking moments and resolving tangentially and radially those forces 
acting on the sliced edges, together with the reversed inertia forces, taken 
over the whole area of the cross-section. Thus 


éo, 
A 
éo, 
J 


where p is the density of the material. The quantities o, and r,, are the 
usual plane polar stress components. 

It is now assumed that the thickness normal to the plane of curvature 
is so small that there is plane stress throughout this thickness and that 
the radial stress component ¢, is negligible in comparison with the other 


(5) 
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components. Thus the stress-strain relations are 
o,=0, 0, == Be, tq = Grr; (6) 


where ¢, and y,, are the circumferential and shearing strains respectively, 
E is Young’s modulus, and @ the shear modulus. The strains are related 
to the non-dimensional displacements 6 and @ in the usual way by 


(7) 


Neutral axis 
Fic. 3. Resultant forces. 


Following Timoshenko’s (10) theory for straight rods it is assumed that 
the displacements are represented approximately by 
6(8,2z,t) = v(8,t)—z(s, t), @(8,z,t) = w(s,t), (8) 

where v and w are the non-dimensional components of displacement at. 
the neutral axis and ¢ is the change in slope of the normal to this line. 

The resultant circumferential force over the cross-section is denoted by 
N, the radial shearing force by Q, and the resultant moment referred to 
the neutral axis is denoted by M, Fig. 3. They are related to the stresses 
by the equations 


1a Oo 6 
' 
Trs + Gs 
R r 
Fic. 2. Forces acting on an elemental slice. | “3 
bs 
os 
( n+QXN ds 
Os 
Q 
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After substituting from equations (6) to (8) and performing the necessary 
integrations, the resultants are expressed in terms of the displacements by 


(10) 
where the quantity Ak? is defined by the expression 


22 dA/(l—z) = — = Ak, (11) 
A A 


after using equation (3). The shear coefficient «, which is introduced into 
the second of equations (10), is determined by taking into account the 
variation of shear stress over the cross-section. It is discussed further in 
section 7 and is similar to the coefficient used by Mindlin (11) in his theory 
of flat plates. 

Similarly, by integrating the a ae of motion (5), it is found that 


The quantities Ak? and Akj which appear in equation (12) are defined by 
the expressions 


— «1—z)dA = AR, ff 2%(1—z)dA = (13) 
a ‘A 


If the neutral axis is a distance eR from the central line then 


Ak? = — [{ (2,—e)dA = eA 


from equation (11), where z,R is a new coordinate measured from the 
central line along the inward drawn normal. Thus the neutral axis is 
a distance Rk? from the central line. Now, we can quickly obtain a 
quantitative appreciation of the quantities k*, k?, and 43, defined in 
equations (11) and (13), by considering their values for slightly curved 
rods when z = z'/R <1. The equation (3) defining the position of the 
neutral axis can then be written 


zdA/(1—z) = 2(l+z)dA = 0, 


: 
= (1+k jw 
= 
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so that from equation (11) we have 


Ak? = — ff ada = ff aa. 


With this definition, the quantity Rk is the radius of gyration about the 
neutral axis; and, since k < 1, this is the same as the radius of gyration 
about the central line. Furthermore, from equation (13) we find that 


2k? (k<1). (14) 


The equations of motion in terms of the non-dimensional displacements 
v, d, w are obtained by substituting for the resultants from equations (10) 
into equations (12). The result of the substitution is that 


c8 


where the abbreviation 


(15) 

(16) 
is used. The matrix notation is used here to show that the coefficients 
are symmetrically disposed about the leading diagonal for, as pointed out 
by Vlasov (12), this is a requirement of the principle of conservation of 
energy. Equations (15) are three simultaneous partial differential equa- 
tions with constant coefficients and they may be solved by standard 
methods. However, it is convenient for the present purpose to derive 


from them a single equation which contains only the single variable w, 
namely, 


E 
= 0 (17) 


és? + 


where the abbreviations 


K = K** = (18) 
are used. 

If the curvature is slight, say 3k <1, then the approximations given 
in equations (14) may be made. Furthermore, the governing differential 


5092 .54 M 
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equation simplifies to 


R\* a 2E\ @ E 
(19) 


when small terms, i.e. terms which are proportionately smaller than (3k)? 
in comparison with unity, are neglected in the constant coefficients. The 
propagation of free harmonic waves is considered by finding solutions to 
equation (19) of the form 
w= Wexp{i(pt—as)}, (20) 

where W is a constant and p is the natural frequency. The wavelength 
L along the neutral axis is given by 

L = 2nR/A, (21) 
and the phase velocity cp of wave propagation along the neutral axis is 
related to the natural frequency p by 

Cp = pR/d. (22 
The characteristic equation of the wave motions, when the curvature is 
slight, is obtained by substituting equations (20) and (22) into equation 
(19) so that we get 
Cp\*{, E 


eG 


4 (Se 4 22\\_ — 0 (ak <1), (23) 
\ eG) 


eG 

which is a cubic equation in (cp/c)*. The roots give the phase velocities 
of three different types of wave-motion and these are radial shear, exten- 
sion of the neutral axis and flexure, in order of decreasing velocities. 

It is shown later that both equations (19) and (23) can be simplified 
further, whilst retaining the same degree of accuracy in the values of the 
constant coefficients, by taking out the terms which control the extensional 
behaviour. 


c 


4. Neutral axis inextensible 
When the neutral axis is assumed to be inextensible equations (7) and 


(8) show that w = éviee. (24) 


Substituting this into equations (10) we find that the resultants Q and M 
can now be written 


er ad 
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and when the circumferential resultant N is eliminated from the equations 
(12) of motion, they become 


(26) 


after substitution from equation (24). These equations of motion can now 
in turn be expressed in terms of v and ¢ by substituting from equations 
(25) so that 


where c is defined in equation (16). Equations (27) are two simultaneous 
partial differential equations and they give rise to the following governing 
equation containing only one dependent variable : 


C R kK ** 2 E\@ 
E\@ a 
28) 
where equation (24) and the abbreviations of equations (18) are used. 


If the curvature is slight, i.e. 3k <1, then the approximations given in 
equations (14) may be made and the governing equation (28) simplifies to 


R\'a,, E 
= 0 (3k <1). (29) 


c 


The characteristic equation of the wave motions is obtained by substi- 
tuting equations (20) and (22) into the partial differential equation (29) so 
that 


2G eG 


(3k <1), (30) 


which is quadratic in (c,/c)?. The smaller root gives the phase velocity 
of the flexural wave. 


| 
_ 
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It is shown later, in the case of slight curvature, that the partial 
differential equation (29), and the characteristic equation (30), are always 
quite adequate for predicting the flexural motions. 


5. Elementary theory 

The accuracy of the elementary theory is often satisfactory when the 
curvature is slight, ie. 3k <1, and the wavelength of motion is long. 
The effects of rotatory inertia and radial shear are ignored and, further- 
more, the extension of the central line is ignored when considering the 
flexural motion and the flexural characteristics are ignored for the exten- 
sional behaviour. In the following derivation, the equations of the 
elementary theory are derived by considering the flexure and extension 
simultaneously by a method similar to that of Waltking (2). 

When the effect of the radial shear is neglected 

$= (31) 

and equations (10) for the resultants in terms of the displacements may 
be written 


(32) 


) (M/R) = — Baw’ 


det)" 
When the effect of rotatory inertia is neglected then the equations (12) 
of motion reduce to 
aN é 1 ew 
after eliminating the resultant Q. As before, these last equations can be 
expressed entirely in terms of the non-dimensional displacements by 
substituting from equations (32) so that 


2 2 


R\?@ ||. 


where ¢ is defined in equation (16). This gives rise to the governing 
differential equation 


/ a 2 


which is identical to that derived by Waltking (2). However, the constant 
coefficients in this equation undoubtedly contain errors which are of the 
same order as (3k)? is in comparison with unity and so the equation can 


= 0, (34) 


“RS 
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be simplified further to give 
(35) 
The characteristic equation of the wave motions is obtained by substi- 
tuting equations (20) and (22) into equation (35) so that 


(1-4 + (SE) =0 (3k <1). (36) 


The roots of this quadratic equation are 


2 


Cp k(A?—1)? 
since 3k <1. Equations (37) and (38) give the phase velocities of the 
extensional and flexural waves respectively and these velocities agree with 
those obtained in the usual way from the elementary theory. The present 
derivation shows again that there is no significant interaction between 
extension and flexure when the curvature is slight and the effects of 
rotatory inertia and radial shear are neglected. 


6. Discussion of the behaviour when the curvature is slight 

As explained earlier, the curvature is considered to be slight when the 
quantity 3k is much less than unity. The quantity Rk is then the radius 
of gyration about the neutral axis (or central line) and there is a simple 
relationship [equation (14)] between the values of k, k,, and k,. Using 
this definition of slight curvature the general governing differential equa- 
tion (17) has already been simplified [see equation (19)] by ignoring small 
terms (i.e. terms which are proportionately smaller than (3k)? in com- 
parison with unity) in the constant coefficients. In the following, this 
process of simplification is extended so as to obtain simple expressions 
for the phase velocities. 

We have just seen that, according to the elementary theory, there is 
no interaction between the extensional and the flexural deformations. It 
is now shown that this is still true even when the additional effects of 
rotatory inertia and radial shear are included, provided that the small 
terms are neglected in the constant coefficients of the governing differential 
equations. 

The effects of rotatory inertia, radial shear, and extension of the neutral 
axis are included in the governing differential equation (19) whereas the 


165 
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neutral axis is considered inextensible for the derivation of equation (29). 


Now, the following relationship exists between these two equations (19) 
and (29): 


(2) 


(3k (39) 
where the leading differential operator is that obtained for the extensional 
motion according to the elementary theory [cf. equation (37)]. It follows 
from this that the extension plays no role in the flexural motion and that 


—0 (3k <1), (29bis) 
is the Timoshenko (4, 5) type equation which governs the flexural motion 
of slightly curved rods. The characteristic equation which corresponds to 
this is 


(3k <1). (30bis) 

Now when the wave number A > 1/3k, i.e. the wavelength L is very 
short and is less than twenty times the radius of gyration about the 
neutral axis [see equation (21)], then the characteristic equation (30) 


reduces to 
E Cp\* E 
— (2) + (2) = ° 


(3k <1 and A > 1/3k). (40) 
This is identical with Timoshenko’s (4, 5) equation for a straight rod 
and it is deduced that for these short wavelengths the flexural motion 
is not influenced by the slight curvature. Furthermore, since the effect 
of the radial shear is most pronounced when the wavelength is short, it 
follows that the value of the shear constant « is adequately determined 
as for a straight rod and so depends only on the shape of the cross-section. 


|_| 
q 
¢ q 
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On the other hand, for smaller values of A the characteristic equation 
yields an expression for the phase velocity: 
Cp\? k2(2—1)? 
The value of A* < 3/k? is calculated from equation (30) by observing, for 


example, that : 
provided that ay/B? < (3k). This means that 


B/x?G)A*(1 +-A?)(A2— 1)? 
which is true provided that A* < 3/k®. Comparing equation (41) with 
equation (38) of the elementary theory we see that the effects of rotatory 
inertia and radial shear are noticeable even for quite small values of the 
wave number A. 
The equation corresponding to equation (41) for the phase velocity cpg 
of flexural waves in a straight rod is 


for small values of A. Now, noting that 
1 4 (1499) +4 


it follows from equations (41) and (42) that 
Cp\? [eps\* 
It is shown below that this simple relationship, which is the same as that 
predicted by the elementary theories for straight and slightly curved rods, 
is true within the limitations of the present theory for all values of the 
wave number A. 

Now the constant coefficients in the above equations were simplified 
by ignoring small terms which are proportionately smaller than (3k)? in 
comparison with unity. In engineering applications the curvature is 
usually considered to be slight for values of k < 0-03 and so these coeffi- 
cients are thus calculable to an accuracy of | or 2 per cent; but this does 
not necessarily mean that the resulting value of the phase velocity has 
the same acceptable accuracy. In fact, this accuracy can only be deter- 
mined with certainty by making comparison with the exact solutions of 
the three-dimensional equations of elasticity—and this is very difficult. 


(3k <1 and At < 3/k*). (41) 


E 


eG 


) (3k <1), 
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However, it is easy to estimate the percentage difference 5, in the value 
of (cp/c) as calculated from equation (43) with that from the original, and 
more general, equation (23). This percentage difference is given by 
28, —(cp/c)*A*k*( E/x?G) 
using Newton’s method, where c,/c is calculated from equation (43) and 


K Kru 


Furthermore, the maximum percentage error 5, arising from the neglect 
of the small terms in the constant coefficients of equation (23) is given by 
100 (cp/e)*| —A+ 2(cp/c)?A’ —3(cp/c)*Atk?( | 
The percentage difference 5, is largest when k = 0-03 and values are given 
in Table 1 for a range of wave numbers A. Since |5,| < |8,| it follows that 
the phase velocity of flexural waves can always be calculated according 
to the simple equation (43) to an accuracy that is well within the limita- 
tions of the present general theory. The value of cpg is calculated from 


Timoshenko’s equation (40), or is determined experimentally for a straight 
rod of the same material and cross-section. 


(44) 


(45) 


TABLE 1 
Values of 5, and 8, for flexural waves (per cent) 


40 


— 2°06 


TABLE 2 
Values of 5, and 8, for extensional waves (per cent) 


A 2 15 


— 0-069 —0°085 — 
0°83 1°38 2°09 


Equation (39) shows that the phase velocity of extensional waves is 
just as predicted by elementary theory, i.e. 
2 2 
1+)? (37 bis) 
c 
The percentage differences 5, and 5, have been calculated from equations 
(44) and (45) and a few values are given in Table 2 for various wave 
numbers A. Since |8,| < |8,| it follows that the accuracy of equation (37) 


A 2 5 10 20 i 30 = ; 
3, —0°059 —0°033 — 0008 0°023 0-021 

| 

q 

i 


ELASTIC WAVES IN A NATURALLY CURVED ROD 169 


is well within the limitations of the present general theory. The dispersion 
caused by the curvature is greatest for the long wavelengths, i.e. for 
small values of A, and it is noted that c, is always greater than the velocity 
ec of long wavelength propagation in an initially straight rod. However, 
although individual waves travel with the phase velocity, the disturbance 


associated with a particular wavelength travels with the group velocity 
(13) Ca where dcp d 


da 
Substituting from equation (37) it is found that 


(2) (47) 


so that, as expected, the disturbance always travels with a smaller velocity 
in a curved rod than in an initially straight rod. 

Finally, as mentioned in the introduction, it is well known in the case 
of straight rods, e.g. Abramson et al. (8), that the restrictive assumptions 
such as are imposed here, and in the elementary theory, upon the deforma- 
tion of the cross-section do not permit a satisfactory investigation of the 
extensional behaviour except when the wavelength is long (i.e. small A) 
compared with the dimensions of the cross-section. In fact, the present 
theory is hardly adequate for a rod of circular cross-section when A > 1/3k. 
Fortunately, however, the effect of the slight curvature is likely to be 
negligible for these shorter wavelengths. 


= (Acp). (46) 


7. Discussion of the behaviour when the curvature is pronounced 

In the case when the curvature is pronounced the simplification of 

equations (14) no longer holds and it is necessary to revert to the more 

general equations. The governing differential equation (17) includes the 

_ effects of rotatory inertia, radial shear, and extension of the neutral axis. 

When a solution is assumed in the form given by equation (20) it is found 
that the characteristic equation is 


c 


where c is defined in equation (16). If the extension of the neutral axis 


q 
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is neglected, then the characteristic equation is 


2G 
lig 1-4-2) exes |—o0, (49 


from the governing differential equation (28). 

It is unlikely that the above equations for pronounced curvature are 
adequate for motions other than flexure where the wavelength is greater 
than the rod dimensions. Furthermore, in the absence of experimental 
results it is necessary to compare with exact theory in order to assess the 
accuracy. Now, there are but few instances where it has been found 
practicable to use exact theory even for the simpler problem of the initially 
straight rod. In the presence of curvature, it is only recently that Gazis (7) 
has published an exact analysis, with numerical results, for the plane 
strain vibrations of a thick-walled circular cylinder. Unfortunately for 
the present purpose, these results are restricted to the low value of the 
wave number A = 2 and so it is possible only to make a very limited 
assessment of the present approximate theory. 

Equations (48) and (49) are valid for plane strain provided that the 


substitution E = Ep (50) 
is made, where Z,, is the increased modulus defined as 
E 
Ep = 1.3 (51) 


with v as Poisson’s ratio. 

The internal radius of the cylinder is denoted by a and the outer radius 
by 6, so that the wall thickness is h = b—a. The radius of the neutral 
axis, calculated from equation (3), is found to be 


= (52) 


(53) 


The quantities /? and k3 are calculated according to equations (13). 

We come now to the determination of the shear constant « which first 
appeared in equations (10). This constant can be adjusted to any suitable 
value without prejudice to the basic physical assumptions, and Mindlin 
(11), in his theory for flat plates, suggested that the value be determined 
by matching the approximate frequency of the axially symmetric trans- 
verse shear motion with that of the corresponding first mode of exact 


and the quantity k* calculated from equation (11) is ' 

b—aj \a 
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theory. The approximate frequency is obtained by multiplying equation 
(48), or (49), by A? and then making the substitution 

Acp = pR (54) 
before putting the wave number A equal to zero to obtain the axially 
symmetric motion. The result of this manipulation is 


Cs 


Elementary theory. 
=—-Equation (58) 
Equation (49 
Equation (48 
Exact theory Gazis) 


Plane strain vibrations, 
c*= B/(1-v? 


o 


Non-dimensional frequency parameter e 


0-2 04 0-6 0-8 1-0 


Fic. 4. Comparison of various theories for the plane strain vibrations of a circular 
cylinder, wave number A = 2. 
from either of equations (48) or (49). In equation (55) the abbreviation 
= Gip 
is made. The exact frequency is determined (7) from the smallest root 
of the transcendental equation 
where J,, Y, are second-order Bessel functions of the first and second kinds 
[see Watson (14)]. The roots of this equation are tabulated by Jahnke and 
Emde (15) and Gazis (7) has plotted frequency curves for the first and 
higher modes. The smallest root is ph/c, = 7 for h/a = 0, which value 
is only slightly higher (less than 10 per cent) when h/a is as large as unity. 
Taking the root ph/c, = 7, it is found by comparison with equation (55), 
that the shear constant is 


2 


The numerical results of the non-dimensional frequency parameter 
phic = Ah) R\(cp/c) as calculated from equations (48) and (49), for Pois- 
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son’s ratio vy = 0-25 and wave number A = 2, are plotted in Fig. 4 for 
values of h/a between 0 and 1. At A/a = 1 the results from equation (48) 
are less than 5 per cent in error from the exact results obtained by Gazis (7). 
Equation (38) of the elementary theory is strictly only applicable under 
certain circumstances when the curvature is slight. Notwithstanding, the 
non-dimensional frequency parameter ph/c given by this equation for the 
cylinder is ph 2A(A2— 1)(h/a)? 
(hja)}? 
and these values are also plotted in Fig. 4. For h/a = 1 this equation is 
18 per cent in error as compared with Gazis’ exact results. 


(58) 
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SUMMARY 

This paper describes small disturbances of an inviscid, electrically conducting 
liquid resting in a circular dish made of insulating material, and subject to a vertical 
magnetic field. The depth of the liquid is assumed to be much smaller than the 
diameter of the dish or than a typical wavelength. The response to an initial 
disturbance is investigated by means of approximate equations which are an ex- 
tension of the linear ‘shallow-water’ theory. The mathematical problem is essentially 
an eigenvalue problem for a fourth-order partial differential equation with mixed 
space-time derivatives. It is found that, as a result of dissipation by Joule heating, 
all the modes of the solution are damped, but the initial rate of decay of the lowest 
asymmetric mode is only about one-sixth that of the axisymmetric modes. Wave 
profiles and current lines for a particular asymmetric motion are shown. 


1. Introduction 

In a recent paper (1) the non-linear and linear ‘shallow-water’ theories, 
which describe long gravity waves on the free surface of an inviscid liquid, 
are extended to the case of an electrically conducting liquid on a horizontal 
bottom in the presence of a vertical magnetic field. The dish holding the 
liquid, and the medium outside it, are taken to be non-conducting. It is 
assumed that the ratio of mean depth to typical wavelength, H,/L = «, 
is small; that the magnetic Reynolds number N (based on the wavelength 
L and on a velocity (gZ)') is O(1); and that the imposed field strength B, 
is such that the magnetic force does not exceed the pressure gradient in 
order of magnitude. These assumptions follow from consideration of 
laboratory experiments with mercury, for which N ~ 4L*! (L in metres) 
and the field strengths are moderate. Approximate non-linear equations 
(appropriate when the wave amplitude and the depth of the liquid layer 
are comparable) and linear equations (appropriate when the amplitude is 
of even smaller order than the depth) are found; their principal feature 
is that the current, like the fluid velocity, is predominantly horizontal 
and independent of the vertical coordinate. That part of the current which 
results from coupling of the fluid velocity with the vertical magnetic field 
(oV \ B, in the usual notation) will not, in general, have a vanishing 

(Quart. Journ. Mech and Applied Math., Vol. XIV, Pt. 2, 1961} 
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normal component at the insulating boundaries of the liquid, and an 
electric field is therefore required (even though the time variations of the 
magnetic field prove to be negligible) in order that the total normal current 
shall vanish. To the natural length scale (L) of the electromagnetic field 
outside, the shallow liquid layer appears as a current sheet. 

The main problem is to determine the fluid velocity and electric field 
strength in the liquid: the external magnetic and electric fields can then 
be found (in principle, at least) from potential theory. 

Although approximate equations are derived in the earlier paper for 
general waves in a shallow, inviscid liquid, they are solved only for certain 
problems of axisymmetric and plane flow, in which symmetry ensures 
that the oV/B, current is in fact parallel to the interfaces, and the 
electric field is negligible to the lowest order. In the present paper we 
take the linear approximation and solve the initial-value problem for 
general waves in a circular dish. The electric field is now not negligible, 
so that all the magnetohydrodynamic variables play a significant part in 
a problem whose boundary conditions can be realized in the laboratory. 
Indeed, the theory now corresponds more fully to the qualitative experi- 
ments of Lehnert (2). Lehnert found that the surface waves of mercury 
in a glass dish 20 em in diameter and about 4 cm in depth ‘disappeared’ 
under the action of a vertical magnetic field. The theory of this paper 
shows that, like the axisymmetric modes considered previously, the asym- 
metric modes are damped, but that the dominant ones decay much less 
rapidly than the axisymmetric modes. The longer life of the asymmetric 
modes is reminiscent of the ‘dynamo problem’ of magnetohydrodynamics 
(Bullard and Gellman (3), Cowling (4)), although there are, of course, 
substantial differences between that problem and the present one (princi- 
pally, that small perturbations about an imposed magnetic field are 
assumed here). 


2. Formulation of the linearized problem 


Let (X,,X,,X,) be cartesian coordinates and let 7 denote the time. 
Consider a flow of conducting liquid bounded by a solid surface X, = 0. 
a free surface X, = H(X,, X,, 7), and the vertical wall Z of a dish. All 
these boundaries are insulating. A vertical magnetic field B, = i, B, is 
imposed; the mean depth of the liquid is H,; and L is a typical wave- 
length. Rationalized m.k.s.q. units are used. The magnetic Reynolds 
number N and a field-strength parameter k are defined by 


N =poligl}~ Ol), k= 
0 


; 
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where « is the magnetic permeability, o is the electrical conductivity of 
the liquid, and p is its density.t 
We seek the leading term in the expansion of each dependent variable 
in powers of « = H,/L. To this end, dimensionless independent variables 
of O(1) are introduced by 
where « is a symbol defined to be 1 in the liquid and 0 outside it. Similarly, 
dimensionless dependent variables of O(1) are introduced by the following 
expressions for the fluid velocity, electric field strength, magnetic field 
strength (or magnetic induction), and current density: 


X = (2.1) 


V = v2, €v5) ) 
E = Bylegl)*(e,, ees) 


B = B,+(e!N B,)(b;, bg, 53) }- (2.2) 
J= (« €)s) 
We also write H = «L(1+), P = epglp, (2.3) 


where P is the fluid pressure. Arguments which justify the orders indicated 
by this scheme have been given previously (1). 

Substituting from (2.1) to (2.3) into the Maxwell equations for curl E 
and cur! B, and into the vertical momentum equation, and retaining only 
terms of the lowest order, one finds that in the liquid 


(2.44, b) 


The remaining equations are also simplified, but before writing them, we 
prepare the way for their linearization by assuming that the initial condi- 
tions introduce a second small parameter, 


The perturbation quantities v,, v2,..., js, 7 in (2.2) and (2.3) are now O(8), 


+ For mercury, with H,/L = 0-1, k = 1-0, and L = 0-1, 0-2, 0-3 (metres), we have, 
respectively, B, = 0-29, 0-24, 0-22 (10* G). 


] 
= 0, —= 0, 
= jy = = jg = (2.5a, b) 
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and, to the lowest order in « and 5, the equations governing the motion 
of the liquid are (explanations follow): 


(2.7) 


Ge, 
OX, Ok, 
+ = 0 (2.10) 
ox, 1 2 ax, 2 


In the horizontal momentum equations (2.8), the last terms represent the 
magnetic force J A B. Since 7 and (e,,e,) are independent of z;, it follows 
that if (v,,v,) are initially independent of x;, as we assume, they remain 
so for all time. Thus all the quantities in (2.7) to (2.10) are independent 
of the vertical coordinate. Equation (2.7) is the usual linearized continuity 
equation. Like (2.4), equation (2.9) states that the time variation of the 
magnetic perturbation is negligible in the Maxwell equation for curl E. 
The last equation, (2.10), follows not from a field equation, but from the 
requirement of zero normal current on X, = 0 and X, = H: it expresses 
the conservation of current within the liquid. 

To these equations we must add the boundary conditions of zero normal 
velocity and zero normal current on Z. If v and 7 denote local coordinates 
normal and tangential to 7, such that (v, 7, 23) form a right-handed system, 


and if suffixes v and 7 denote the corresponding components of a vector, 
these conditions become 


v,=0 and e+v,=0 (2.11 a,b) 


and, by (2.8), the second may be written 
dv, on _ 


In what follows, v and e denote the vectors (v,,v,) and (e,,¢,), and V 
the vector operator (@/@2,,0/@x,). Taking the curl of (2.8), and using 
(2.10), we obtain 2 
AV) = 0, 


so that if v is irrotational initially, as we assume, it remains so. Then 
(2.10) reduces to 


V.e=0. 


3 
: 
4 
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We introduce a velocity potential ¢, and an electric vector potential i, x, 
such that by ax 


and the governing equations become 


(2.12) 
= —Vy—2kV(x+4), (2.13) 
V2y = 0. (2.14) 
Taking the divergence of (2.13), we obtain upon substituting for V*¢ from 
2.15 
= 0, (2.15) 
or operating with @/ét, 
é 
(v 0. (2.16) 


In view of the boundary conditions, the problem cannot be formulated in 
terms of », and we take (2.16) as our governing differential equation. 
The boundary conditions (2.11) become, when (2.11¢) is differentiated 
with respect to time, 


As initial conditions, we take 
= F*(x,,2,) and = G*(z,,%,) att=0, (2.18a,b) 


where 0F*/év = @G*/év = 0 on Y. The arbitrary constants in the poten- 
tial functions F*, G*, y, and 7 at t = 0 are fixed by the requirement of 
zero mean value: 
ff F* dz,dz,=0, ete. 
The initial values of y and 7» are then related to F* and G* by 
= 0, -= on 


= + F*) 
(If v = 0 initially, then F* = y = 0 and n = —G* at t = 0.) 

The problem presented by (2.16) to (2.18) bears some resemblance to 
classical vibration problems, such as those for plates, but differs signi- 
ficantly from these in the appearance of mixed space-time derivatives in 
the differential equation, and of a time derivative in the second boundary 


5092 .54 N 
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condition. In particular, we cannot expect the eigenfunctions of the solu- 
tion to be wholly orthogonal. 

It is not difficult to show that any solution of (2.16) may be decomposed 
into such that 


ale =0 and V%,=0.  (2.19a,b) 


3. Eigenvalues for the case of a circular dish 

Let (r,@) be polar coordinates corresponding to (z,,x,), and consider a 
circular dish of radius L; then Z becomes the circle r = 1. If 6 = e*%9(r, t) 
(ns = 0, 1, 2,...), we have 

2 12 # 
= € Die, say. (3.1) 
Frequent use will be made of Bessel’s equation, 
(D3+-a*)J,(ar) = 0 or (D2—6*)I,(br) = 0, 
where J,, and J, are, respectively, the Bessel function and modified Bessel 
function of the first kind and of order n. 

For the axisymmetric modes (n = 0) the electric field vanishes, since 
it is the gradient of a harmonic, axisymmetric potential and must be free 
of singularities in r < 1. Hence we merely have to solve (2.19a) with the 
boundary condition (2.17a), and the appropriate eigensolutions are 


$ = (3.2) 
where A,,,, is the mth positive zerot of J‘,, and w,,, = (A2,,—k*)'. The 
Wom are all real for k < Ag, = 3-83. 

For the asymmetric modes (n > 1) the elementary solutions of (2.16) 

which are finite at the origin are 

(qr) and (3.3) 

where p is a complex number, and gq? = p?+2kp. For the sake of definite- 

ness, we draw a cut from p = —2k to p = 0 and let ¢g denote that branch 

which is positive on the positive real p-axis: later the choice of branch 
will be seen to be immaterial. The solutions 

= I, = k,n), say, (3.4) 

satisfy the first boundary condition, (2.17a); the second, (2.17b), then 

requires f(p; k,n) = 2kn = 0. (3.5) 


The eigenvalues p are the non-trivial roots of this equation; by non-trivial 
roots we mean those not at ¢g = 0 (p = 0 or —2k). The form of f shows 


+ This will be written A,,, only where there is a possibility of confusion: e.g. when 
m = 


t 
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that, for k real, the complex zeros occur in conjugate pairs. The corre- 
sponding eigenfunctions ¢ also occur in complex conjugate pairs. 

In some contexts it is useful to consider, instead of f, the function 
f* =f, where B is 1 if n is odd, and 0 if n is even: f* is an analytic 
function of p and k for all finite values of these variables, and f and f* 
have the same non-trivial zeros. 

We first establish two general results about the zeros of f, which are 
stated as lemmas. 

Lemma 1. Let C, be a circle |p| = R, in the p-plane, such that 

Aneto < (a > 0, = 1,2,3....). 
For k and n bounded, there exists a finite number A = A(a,k,n) such that 
for all R, > A the function f *(p; k,n) has 2s non-trivial and n+-B+-1 trivial 
zeros within C,, a zero of order | being counted | times. 

To prove this, we first note that near g = 0 

~ Ofp.[p(p+ 2k) }"*?}, 
and f* ~ 
so that a circle C, for which R, > 2k contains n+8-+-1 trivial zeros of 
either function. Next, we observe that a circle |p| = A,, > k passes 
through the points p = —k+i(A2,—k*)', which are the sth pair of non- 
trivial zeros of pq'+#I’,(q). A circle C, for which R, > 2k therefore contains 
2s+-n+8+-1 zeros of this function. 

Now, by Rouché’s theorem,t pq'+*J',(q) and f* have the same number 
of zeros within any contour on which 


> |2kng®l,(q) |, 
or \pq| > |. 
But with p bounded away from the zeros of /',(q), the right-hand side is 
bounded, and this inequality is true if R, is taken sufficiently large. This 
completes the proof. 

The lemma shows that the eigenvalues form an enumerable set, and 
that those having a modulus greater than some finite value lie near the 
circles |p| = A,,. 

Lemma 2. Fork > 0, and with k,n and |p| finite, the complex eigenvalues 
lie in —k < rep < 0, and the real eigenvalues (if any) lie on —2k < p < 0. 

The proof is based on a standard method of eigenvalue problems. The 
eigenfunctions #, defined in (3.4), satisfy the equations 

= 9, (3.6) 


+ If ¢(z) and ¢(2) are analytic within and on a closed contour C, and |¢(z)| > |¥(z)| on C, 
then ¢(z) and ¢(z)+(z) have the same number of zeros within C. 


a 
i 
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and (Dy) = —2kpmy atr=1, (3.7 a, b) 


the last of which follows from (3.4) and (3.5). Multiplying (3.6) by J, the 
complex conjugate of ¥, integrating from 0 to 1, and integrating by parts, 
one obtains 


1 
0 


1 
| dr=0. (3.8) 
0 


The integrated terms vanish at r = 0, and those at r = | are simplified 


by means of (3.7). Let 
1 


| dr = A, > 0, 
0 
1 


1 


i 
2 n2 
0 

Here we have been able to write > rather than > because it is easily 
verified from (3.4) that the relevant terms do not vanish. Equation (3.8) 


now becomes A,+2kpA,+p*A, = 0, 
A A? A,\} 
P A,+\ A 


Since A, > A, > 0 and A, > 0, this proves the lemma. 

Accordingly, the asymmetric modes are damped, but in the absence of 
real eigenvalues (that is, in some neighbourhood of k = 0) they are not 
damped as strongly as the axisymmetric ones, (3.2). 

It is now clear that, as & is increased from zero, the eigenvalues move 
from the points p = -+-iA,, into the half-plane rep < 0. Each eigenvalue 
p is a single-valued function of k by definition, and its derivative 


is uniquely defined in a complex k-plane,t except possibly at those values 


+ We shall continue to treat k as real and positive except where the contrary is explicitly 
stated. 
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of p and k for which 
f*=0 and @*/ap=0, (3.9) 
that is, at non-trivial double zeros of f* or f. Away from such points, 
the eigenvalues p are therefore analytic functions of k. In view of Lemma 
1, a non-trivial double zero of f* in the p-plane can occur only by the 
confluence, as k increases, of two of the eigenvalues. Moreover (for k real 
and positive), such a confluence of an eigenvalue and its complex conjugate 
must occur, on the negative p-axis, before real eigenvalues are possible. 
The roots of (3.9) are therefore of some interest. Eliminating J*(q) from 
éf/ép by Bessel’s equation, and J',(q) by the equation f = 0, one obtains 
At a non-trivial zero of f, I,(q) cannot vanish.+ and the conditions 
f=0 and =0 (3.11) 
are therefore equivalent to (3.9). 

For a given n, let k, be the complex &-root of (3.11) having the smallest 
modulus; let k, be the smallest real, positive k-root; and let k, be the 
smallest real, positive k-root for which the p-root is real. (It is, of course, 
possible that two, or even all three, of these k-values are equal; and, indeed, 
it is probable that k, = k,, for k, merely corresponds to complex con- 
jugate q-roots on req = 0, img < k, whereas k, + k, implies a confluence 
in both the qg- and p-planes away from lines of obvious importance.) Then 
for 0 < k < |k,| ali the eigenvalues may be represented by convergent 
power series in k; for 0 < k < k, no double eigenvalues exist and all the 
eigenvalues are analytic functions of k; and for 0 < k < k, no real eigen- 
values exist. Values of k, and of the corresponding p have been computed 
for n = 1, 2, 3 and are shown in Table 1: details of the calculation are 
given in Appendix A. For n = 1 it has also been confirmed that k, = ky, 
and it seems reasonable to assume this for all n. 


TABLE 1 


hy 

1 4°93 —451 
2 6°07 —5'51 
3 7°23 —663 


For k < k,, we may now label the eigenvalues p,, (s = ..., —2,—1, 
1,2,...), where p,, is that root of (3.5) which > +iA,,,, as k-> 0 (+ for 


+ For this would imply either pg = 0 (a trivial zero) or I/,(¢) = 0, which is impossible 
at J,(¢) = 0. 


: 
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8 = 0): the corresponding eigenfunctions are labelled %,,,. To compute the 
Png» We derive a power series in k. With A,,, = A (s > 0) and the substitu- 


Png = At+yk+iy, + Ro, (3.12 a) 
1 


the equation (3.5) is expanded about qg = iA: all the derivatives of J,,(q) 
may be expressed, by means of Bessel’s equation and J',(iA) = 0, as 
products of J, (iA) and algebraic functions of A. Writing 


2n 
¢, 


(3.13) 


we obtain 


Ye = Fc? — hnc*) (3.14) 


= 
These coefficients are shown in Table 2 for values of n and s up to 3, and 
it is clear that the series (3.12a) has every appearance of convergence for 
0<k < 2, say. This is probably sufficient for practical purposes, since 
k > 2 corresponds to rather strong magnetic fields. The series, truncated 
after y,k*, has been tested in (3.5) for n = 1, s = 1 (the worst case), and 
k = 1-0; the result was |f| = 0-0005. 


= —I+¢ | 


TABLE 2 


18412 8°5363 
—0°1632 —0°9722 
0°0035 —o'0617 
o"ooo! 


3°0542 9°9695 
—0°2493 —o'9581 
00058 


4°2012 11°3459 
—0°3064 —0°9499 
0°0068 —o'0481 
—0°0012 


Certain simple properties of y, can be derived analytically, and confirm 
the trends shown in Table 2. Since 
Mie > > 
[(5) 486], y, decreases (algebraically) as s increases with n fixed, and is 
always in the range —1 < y, < 0, as would be expected from Lemma 2. 


| = 
| 
n=1 A 
& v1 
Ys 
Ys 
n= 2 A 
Ys 
n= 3 A 
3 "1 
Ys 
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For n = 4, 8 = 1, we have y, = —0-3488, and for n > 4 [(5) 487] 
2, >n?+3n, sothat —1 <y, < —}. 
Hence the value of y, for n = 1, s = 1 is the largest. 

While the series (3.12) have been derived on the basis of k > 0,n ~ O(1), 
and s ~ O(1), they are also useful for k ~ O(1), n> and/or s >. 
This might be expected, since p and q are large in this case, so that the 
second term of f still dominates. From the properties of Bessel functions 
of large argument (5), and of large order (6), we find that 

n 
Ans ~ O(n+8), c~ (3.15) 
and it follows that the series (3.12) contain the first few terms of asymp- 
totic series with respect to n or s (as well as terms comparable with, or 
smaller than, the remainders). For (3.15) shows that the terms in (3.12a, b), 
up to those in k*, successively decrease in order of magnitude for k ~ O(1), 
nm and/or s-» 00. Moreover, if the derivation of the series is re-examined, 
and the remainders are carefully estimated at every stage, one finds that 
Ry ~ for s < O(n), 
~ O(k'ns-*) for s > O(n), 
which is always smaller than the dominant part of the k* term in (3.12b), 


and that i 
Ry ~ — + 


4. Solution of the initial-value problem 

We shall proceed formally, and establish sufficient conditions for the 
rigorous justification of our solution a posteriori. 

The initial-value functions F*(x,,2,) = F(r,@) and G*(x,,2,) = G(r, 6) 
are first expanded in a Fourier series with respect to @ and a Dini series 
(5) with respect to r. Thus 


F(r,0) = re > nm (Aum?) 


where 


Fum = 


and 4, is 1 for n = 0, and 0 otherwise; G is expanded similarly. 

The orthogonal functions J,,(A,,,, 7) are not eigenfunctions of our problem 
for n > 1, but it will appear that they form a convenient stepping-stone 
from the arbitrary functions F and G to the non-orthogonal functions 
,,(r). It is also useful to bear in mind the ‘fully damped solution’ ¢™, 


} 


n®) (Anm) } J Fir, r)r 


4 
(4.1) 4 
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which results from neglecting the electric field. This satisfies 

= 0, 0 om? = I, (4.2) 

and the same initial conditions as the required solution: it is 


sin w,,,, t)+- 


+ Gyr 810 t}Tn(Anmt), (4-3) 
where w,,,, = (A2,,—*)!; (ef. (3.2)). 
Returning to the real problem, we pose a solution of form 


and — Gy Jy Quy?) att = 0, (4.4) 


and introduce the Laplace transform 


= p | e-Pig(r, t)dt. 
0 


Here p is a complex variable which will be identified presently with that 
of section 3; as before, we write p?+-2kp = q?. The transformed problem 


for each 9, (with the suffixes on ©, Anm> Fam, and G,,, temporarily 
omitted) is then 


(Di—P)Dne = MPF (Ar), (4.5) 


atr= (4.6 a) 
or 


and (D2—p*)6 = —(p*F+pG@)J,(A) atr=1, forn > 1. (4.6 b) 
A particular integral of (4.5) is given by 
go — 


and this also satisfies (4.6 a): in fact it is the transform of the nm-component 


of the fully damped solution. For n = 0, @ suffices; but for n > 1, we 
must add the complementary function 


J, (Ar), (4.7) 


2kJ, —pG){nl, (qr) (4.8 a) 
__ (AYP F —pG)p(r, p) (4.8b) 


(q*+-A*)f(p) 


24 

n=0 m=1 
n=0 m=1 
where 


in order that both the boundary conditions (4.6, b) be satisfied. Then 


Pp 


+ t+ (49) 


where ¥ is the Bromwich path C—ioo to C+ ico (C > 0). We observe 
that ¥ and f are both even or both odd functions of ¢g, so that a circuit 
in the p-plane about a branch point of g leaves the integrand of (4.9) 
unchanged, and this integrand is a single-valued function in an uncut 
p-plane. Also, the points g = +iA,,, are not poles of the integrand, 
because the expression in { } vanishes there. Nor are the points p = 0 
and p = 2k poles, for there 4 ~ Of{pf(p)}. Thus the only singularities of 
the integrand are poles at the non-trivial zeros of f, that is, at the eigen- 
values p = p,,: these are simple poles except for certain points at parti- 
cular values of k > ky. 

The integral around a large semicircle in rep < 0 vanishes for t > 0. 
Evaluating the residues at the poles, using (3.10) for f’(p,,), and summing 
over n and m, we find that 


n=1 m=1 8=—@ 


A Pns (Anm) B Pre 2kJ),(Anm) 
(Gne Gast Past n(Ins)( Gas + rm) 
(4.104) 


Pras 2)p—2k*n? (p Pns)s (4.10e) 
and ¢” is the axisymmetric (n = 0) part of the fully damped solution. 
Sufficient conditions for the validity of the solution (4.10) are given in 
Appendix B.t 
Because of the nature of the problem, the expressions (4.10) are some- 
what involved: we now try to clarify their meaning. We shall assume 
that k < 2 (say), so that (3.12a) and Table 2 give acceptable values of 
the p,,,, and, in particular, y, in Table 2 shows the rate of decay of the 
various terms. Consider an s-series, with n and m fixed: we may regard 
this as a ‘mode’ of the motion in the sense that it expresses the response 
to a single component of the initial disturbance. Two terms of the s-series 


+ A uniqueness proof of the usual type follows directly from the energy equation (5.2), 
integrated with respect to time. 
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are of special importance: that with s = 1, which dominates at large 
times, and that with s = m, which is the only non-vanishing term for 
k—+>0 (Gy. iA,,) and is still strongly dominant for small times when 
k = 1-0. For example, with F = 0, r = 1, and t = 0, (4.10) and (4.4b) 
yield 


Prone 
1 (for all n, m), (4.11) 
and with k = 1-0, n = 1, and m = 1, 2, or 3, the difference between the 
8 = m term and the total is less than 2 per cent. (The error after 3 terms 
of the series is less than 0-1 per cent.) 

Alternatively, we may prefer to visualize another type of ‘mode’: that 
given by an m-series with n and s fixed. Then the m-series merely repre- 
sents a constant amplitude and phase shift for a single eigensolution, and 
the mode is characterized by one rate of decay and one frequency. Of 
these modes, that with n = 1 and s = 1, which varies in time approxi- 
mately as exp{(—0-16k-+-11-84)t}, decays most slowly and has the smallest 
frequency. However, this viewpoint also has disadvantages: a mode of 
the second type is not directly comparable with its fully damped counter- 
part or its counterpart for k = 0, since its shape depends on the complete 
initial conditions in a more complicated way; and it has no convenient 
‘small-times’ representation like that found below for modes of the first 
type. 

Neither type of mode has the usual property of preserving a fixed spatial 
shape which is merely scaled up and down with varying time, but at a 
fixed point the second type of mode does vary exactly like a damped 
oscillation. 

Since the velocity potential is known, the surface elevation 7» and the 
electric vector potential y are readily found from 


t 
n(r, = n(r, 0, 0)— dt, 
0 


and 


If Hnm and yx,,,,, are defined by equations like (4.4a), we have 


Tom = 24} 
x (4.12) 


‘am = — (Gna?) > 1) 


s=-@ 


and Xam = —™™ > 1). (4.13) 
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In the remainder of this paper we shall set F = 0 to avoid unduly long 
expressions: the initial values now correspond to zero fluid velocity, zero 
electric field, and to a surface elevation » = —G. 

Another representation of the solution is possible: this is useful for small 
values of the time and is of some physical interest. We split each mode 
of the first type, ¢,,,, into three parts: (i) the fully damped part, repre- 
sented by the J,(A,,,,7) term of (4.9), and corresponding physically to a 
wave in an unbounded fluid (in which case no electric field would be 
required); (ii) a part represented by the J,(qr) term, corresponding to a 
wave which is associated with the »-field and which propagates from the 
wall with dimensionless velocity one (and dimensional velocity (gH,)*); 
and (iii) a part represented by the r” term, corresponding to a wave which 
is associated with the electric field and which propagates from the wall 
with infinite velocity. Evaluating the inversion integrals of (ii) and (iii) 
by expanding for p + 0 (that is for t-> 1—r and for ¢ > 0, respectively), 
one obtains for n > 1 


Pam SiN Want In(Anm?) (t+r—1) x 


Gum nm 


where U is the Heaviside unit function. Such series are not asymptotic 
uniformly with respect to n, m, and r, but they are useful for individual 
modes 9,,,, When n, m, and r-! are bounded. 

The corresponding expressions for 7 and y are 


— 


xU(t+r—1) (4.15) 


and 


5. Results and discussion 

The principal result of this investigation is that, whereas the axisym- 
metric modes of the solution vary in time like 


exp{[ —k*)*}t}, where Ay, = 3-83, 


i 
| 
| 
| 
| 
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the dominant (n = 1, s = 1) asymmetric mode varies approximately like 
This result may be partly explained in terms of the energy balance,t 
which is also closely related to the proof of Lemma 2. 
Using (2.4) and (2.5) we write for the (dimensionless) current j and 
Poynting vector s in the liquid 
08. 
(EAB, makes no contribution because V(E A B,) = B,(V AE) = 0 to 


our approximation, and the flux of (s,,8,) is across small surfaces in the 
liquid and negligible.) Since (2.12) and (2.13) yield 


ov 
= x), 


we obtain for a volume of fluid, bounded by a vertical cylinder C of cross- 
section A, 


d 
nt dr = at $v?) dA+ 2k If 44, (5.1) 
Cc A “A 


where v and 7 are coordinates normal and tangential to C. This energy 
equation states that the rate of doing work by pressures on C balances 
the rate of increase within C of potential and kinetic energy, the dissipa- 
tion by Joule heating, and the emission of electromagnetic energy from 
the upper and lower surfaces of the liquid. However, this last part of 
the energy is merely communicated to other parts of the liquid, and not 
to the external field: this can be seen in at least two ways. 

The external field is quasi-steady to the lowest order (curle = 0) and 
current free (curl b = 0), so that 

divs = div(e Ab) = 0, 

and no part of the external field absorbs energy. The rate of change of 
magnetic energy, 0(B*/2u)/é7', is of higher order by a factor O(e) than 
the energy changes represented in (5.1). Again, if we choose C in (5.1) to 
be the dish 9, of cross-section .o/, the s,-term is (since V?y = 0) 


JJ = ($+x) ar, 


+ In the previous paper (1) the energy equation was written correctly, but the significance 
of s, was not interpreted properly. 
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Now by the boundary condition of zero normal current, (2.11 b), (¢+ x) 
is independent of tr on Z, so that this expression is proportional to 


| 
9 


Since @¢/@v = 0 on &, the energy equation for the whole dish becomes 


| = —2k | pda. (5.2) 


Thus the total mechanical energy is continuously dissipated by Joule 
heating. Now consider (i) an axisymmetric mode, for which the current 
lines are circles concentric with J, (ii) a low asymmetric mode of the fully 
damped motion (which may be visualized as that in a liquid unbounded 
laterally), and (iii) the corresponding asymmetric mode of the actual 
motion. The average current densities of (i) and (ii) are of comparable 
magnitude (the damping is the same for the two cases), but the average 
current density of (iii), at early times, is substantially smaller. This is 
because in the unrestricted case (ii) j is normal to v, and the current lines 
cross Y normally, whereas in case (iii) the insulating wall has been placed 
across these ‘free’ current paths. In view of (5.2), these lower current 
densities (at early times) lead to a slower decay of the actual asymmetric 
motion. (For the high modes, the current in any case circulates within 
a large number of small cells, and the wall has a smaller overall effect.) 

It remains to consider the distributions in space of the various fields, 
and the phase differences between them: these questions are perhaps best 
treated by means of an example. We take k = 1-0 and the simple case 


8 cos 7) 


F = 0, G(r,8) = —n(r,0,0) = — FO)? 
M1 


that is, we study the mode n = 1, m = 1, with G,, = —8/J,(A,,). 

The figures show the surface elevation, 7, and the current function, 
$+ x, as representative of the hydrodynamic and electromagnetic distur- 
bances, respectively. The computations were based on the s-series of (4.10 b), 
(4.12), and (4.13). It was sufficient to use (3.12) for the p,, and to retain 
only terms up to s = 3: checks were provided by the values at t = 0, the 
maximum error (at r = 0-3, 8 = 0) in n/8 being 0-005, and that in (¢+ x)/ 
being 0-001. For Fig. 2 both the s-series (4.12) and the ‘small-times’ 
approximation (4.15) were used: the results are indistinguishable to the 
scale of the figure. 


§, 
| 
‘ 
id 
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Fig. 1 shows the variation with time of the surface elevation at the wall 
of the dish; the fully damped solution and the exponential decay rate of 
the first (s = 1) component are also shown. This component dominates 
to such an extent that the departure from a damped, purely sinusoidal 


1:0 


Fic. 1 The surface elevation at r = 1. present solution, -~-—— fully damped 
solution, exp(re ¢). 


i 


0 
9 0-2 : 0°6 0-8 i 1-0 


Fic. 2. Surface wave profiles at small times. ——-— present 
solution, -—-—-— fully damped solution. 


motion is very slight. Figs. 2 and 3 show wave profiles at small times and 
at times which correspond approximately to quarter-periods of the surface 
motion at the wall. The difference between the complete solution and the 
fully damped one (Fig. 2) is the incoming wave in r > 1—t which was 
mentioned in section 4. It is clear that the complete wave profile changes 
its shape with time, the elevation near r = 0-5 lagging somewhat behind 


| 
cos 
0 
3 4 
-0°5| 
4 “10 
1-0 = 
0-2 
cos 
0:4 
0-2 
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T 


t=0 


. Surface wave profiles at t = 0, 0-8,..., 3-2. 


fame 


Fic. 4. Current lines ¢+ = constant. The increments of (¢-+-x)/5 are 0-016. 
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that at r = 1, but this is not a dominant effect and the profile shape after 
one period resembles the original one. At the wall there is a radial pressure 
gradient which is balanced by the magnetic force, and at the wall this 
force stems only from the azimuthal component of the electric field: 


jo=% = 


Thus Fig. 3 suggests that the current density at the wall leads the surface 
elevation there by a phase angle not very different from }7. 


T ' 


i 


Fic. 5. The current function $+ y at r= 0-5, 8 = 0. 


Fig. 4 shows current lines at t = 0-8 and t = 2-4; the former is near the 
time of maximum current density throughout the dish. Near the centre 
of the dish, and near the ,- and z,-axes, the current is essentially normal 
to the fluid velocity, as it would be in an unbounded fluid according to 
our approximation, but of course both fields are tangential to the wall. 
It appears from further calculations that the shape of the current lines 
varies little with time, and that the points of zero current, where |¢+ x| 
has maxima, remain between r = 0-5 and r = 0-6 on 0 = 0 and z. Thus 
o+yx at r = 0-5, 6 = 0 is a measure of the total current circulating in 
each half of the dish: it is plotted against time in Fig. 5. Broadly speaking, 
we may say that the current as a whole leads the surface elevation by a 
phase angle of approximately }7. 
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APPENDIX A 
COMPUTATION OF k, AND &, 


A change of notation is convenient so that certain symbols have different meanings 
in this appendix and in the main text, but no confusion should result. With 


(3.11) becomes s*—xs* + 8*+ = 0, (Al) 
and 2nJ,(nt)—stJ) (nt) = 0. (A2) 


« is assumed real and positive, and (A 1) then has the following s-roots: s;, real and 
such that 0 < sy < «; 8, real and negative; and sy, sy, a complex conjugate pair 
which lies in 0 < res < « for x? > 4(N+1). Hence s,, sy, and syy are possible 
double roots of f. 

To find the k, of the text, that is, to find the smallest «x (= «;, say) such that 
8 satisfies (A 2), we write x = vs, and (A 1) yields 


1<v < HN+(N*+8)}} < 1-78. (A3) 


Taking values of v which descend from }{N+(N*+8)#}, we conipute s; and hence 
t;, and test it in (A 2), which may be written 


= 0. (A4) 


The largest value of v to satisfy (A 4) corresponds to the smallest «x, that is, to xs. 
The question is now whether for 0 < x < Ks, 8yy satisfies (A2); if not, then 
k, = ks, as conjectured in the text. 8, is computed from (A 1) by 


+ 
= 


v < HN—(N*+8)}}, 0 < way < Ks, 


and, now that s; and sy are known functions of x, the complex roots of (A 1) are 
easily found. These sq, and ty are then tested in (A2). The complex Bessel 
functions can be computed fairly quickly from Taylor series for n = 1, when 
im(nfyyy) < 1-5, but the computation would become laborious for n > 2, when 
im(ntyzy) is larger, and has therefore not been attempted for these cases. For n = 1 
and 0 < « < x, it was found that sy did not come close to satisfying (A 2). 


APPENDIX B 


THE CONVERGENCE OF CERTAIN SERIES 
DERIVED FROM (4.10a) 


To justify the solution (4.10a) we split it into two parts, as in (2.19): dy consists 
of ¢\? and of the terms involving [,(¢,,7) in the triple series, and ¢yq consists of 
the terms involving r*. Then if two term-by-term differentiations with V or @/ét 
can be justified for each part, it follows that ¢y satisfies (2.19a) and ¢y (2.19b), 
and that the boundary conditions and initial conditions are satisfied. Sufficient 
conditions are given by the following. 


5002 .54 


4 
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3. If,inr < 1, F*(x,,2,) = F(r,@) and G*(x,,x2,) = G(r, 0) have continuous 
derivatives with respect to (x,,%_) up to the second and first orders, respectively, those of 
the next order being of bounded variation, then the series obtained from (4.10a) by up 
to two termwise operations with V or 2/ét converge absolutely and uniformly in r < 1 
and t > 0. 


In the following, a convergence proof is outlined for the terms involving G,,, in 
the triple series; the proof for the terms involving F,,, is entirely similar. Since we 
are establishing absolute convergence, it is sufficient to take s > 0 and to consider 
orders for n + co and/or m > © and/or s > 0; we shall so express our estimates 
that they are uniformly valid for all n > 1, m > 1, and s > 1, and for all r < 1 
and t > 0. 

From (5) and (6) we find that 


Anm ~~ O(n+m), ~ Of{mi(n + m)t}, (Blab) 
r)~ O(n-* Ofm-*(n m)~*}. (B2ab) 
(i) We first have to estimate the G,,,, which are defined by 
no 
G,(r) = G(r, (B3) 
1 
2)2 
Gam = G,(r)J,(Ar)r dr, (B 4) 
0 


where A = A,. In (B4) we integrate by parts twice, apply recurrence formulae for 
and J,,,(A), and use G)(1) = 0, to obtain 


1 1 
G,(r)J,(Aryr dr = | (B5) 


where p = }r* and g, = r-"G,. It may be shown that, if ¢-*h(t) has bounded 
variation in [0,1], then 


{ h(t)J,(At) dt ~ O(vta-4). (B6) 
This is an extension for v + © of the result on p. 595 of (5). 

To establish the bounded variation of r~*(d*g,/dp*)r"**, particularly near r = 0, 
we observe that the restrictions on G(r,0?) = G*(x,,2,) permit us to write 


G(r,0) = G*(0, 0)+r{cos 6 GF (0, 0)+sin 6 GF (0, 0)}+ (B7) 
0 


where suffixes denote partial derivatives. The last term is O(r*?), of bounded varia- 
tion with respect to #, and has two derivatives with respect to r, of which the 
second has bounded variation. It follows from (B7) and (B 3) that r-#(d*g,,/dp*)r"** 
has bounded variation in [0,1], and is O(n~") for n + oo, 

Using (B 6) and G,(1) ~ O(n-*) in (B 5), we then have 


G,(r)J,(Aryr dr ~ Ofm-*(n + 
6 
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and the second term on the right dominates. Hence in (B 4) 
Gam ~ O{n-tm-t(n (n > 1) 
~ O(m-4) (n = 0) 
(ii) We now discuss the absolute convergence of derivatives of 


(B8) 


n=l m=1 


1 
= — ns - 
Pnmsl nm n(Anm) P,, Gas + 
Ins 

and P,, is the quartic defined in (4.10¢e). Since asymptotic formulae for the p,,, 
and q,, (with k ~ O(1), n and/or s > 0) are given by (3.12), terms like 9,1 are 
estimated without great difficulty, but care must be taken with P,, and g3,+Aim 
in which cancellations between the dominant terms can occur (for p,, ~ im and 
8 = m, respectively). 

We shall consider D*9,,,,1, where D® denotes D2 or D® Cte. 
have orders equal to, or smaller than, those shown. It is necessary to divide the 
‘space’ n > 1, m > I, # > 1 into two parts: the ‘plane’ s = m, and the domain 
jm—s| > 1. One finds that on s = m, 

~ Ofn-tm-t(n +m)-4} 
< O(n-tm-4), (B9) 
which ensures absolute convergence: while in |m—s| > It 


Pnmll = Pnmsl 


~ + m)-4 |m*—s*|-(n + 8)} 
< (m2 — 8? |-1), (B10) 
The absolute convergence of these terms is amply proved by reference to the 


integrals dmds st dmds 
A 


where A is the domain m > 1, ¢ > 1, |m—s| > 1: these integrals converge. 

(iii) It remains to estimate the axisymmetric (n = 0) terms: at the same time 
we consider the entire fully damped solution (4.3). From (B 1), (B 2), and (B8) one 
finds that O(m-4), 
and ~ (n > 1) 

< O(n-im-4) 
(cf. ON 8 = M). 


+ The estimate (q2,+23,,)"? ~ O(|m*—s*|-1) is conservative for m ~.o(n), 8 ~ o(n). 
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AN APPLICATION OF MANGLER’S TRANSFORMATION 
By N. RILEY (Department of Mathematics, University of Manchester) 
[Received 18 August 1960] 


SUMMARY 
Using Mangler’s transformation, wall jets on bodies of revolution are investigated. 
The analogous two-dimensional flow occurring is found to be the two-dimensional 


wall jet, independent of the particular body of revolution. Two examples of such 
wall jets are given. 


1. Introduction 

MANGLER (1) has shown that in the absence of swirl there is a relationship 
between axially symmetrical and two-dimensional boundary layers which 
reduces the calculation of the laminar boundary layer on a body of revolu- 
tion to that on a cylindrical body. For a given body of revolution the 
analogous two-dimensional flow can be determined. The two-dimensional 
analogue of an axially symmetrical boundary layer depends both on the 
pressure distribution and on the shape of the body of revolution; thus an 
axially symmetrical stagnation point flow corresponds to symmetric two- 
dimensional flow past a wedge of angle $7. 

If we consider the problem of a circular jet impinging on an arbitrary 
body of revolution along its axis, Mangler’s transformation shows that the 
problem can be reduced to that of a plane jet striking a plane surface 
normally and spreading out over it. Glauert (2) first investigated such a 
two-dimensional wall jet and also examined in detail the radial wall jet 
formed when a circular jet impinges normally on a plane surface and 
spreads out radially over it. Using Mangler’s transformation the present 
paper investigates wall jete on bodies of revolution, Glauert’s radial wall 
jet occurring as a special case. In particular it is shown that for these 
wall jets the flux of exterior momentum flux is constant and the additional 
parameter introduced by Glauert (3), using techniques due to Watson (4), 
is recovered. It may be noted that if the body is convex to the fluid then 
the flow may separate before the final similarity form is attained, in which 
case the results derived below no longer hold. 

Mangler’s theorem also holds for compressible boundary layers and so, 
provided that temperature differences are small enough for the effects of 
gravity to be negligible, compressible wall jets on bodies of revolution 
may be studied using the results obtained by Riley (5). 


+ Now at the Department of Mathematics, Durham Colleges in the University of Durham, 
Durham. 
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2. Application to wall jets 
The boundary layer equations for si ve two-dimensional steady 


where # and g are measured along and normal to the boundary respectively, 

a and @ are the velocities in these directions, and U,(2) is the external flow. 

The corresponding equations for axially symmetrical boundary layer flow 

without swirl are U atu 


(4) 


(ru) + 

where again x and y are measured ti and normal to the wall, « and v 
are the velocities in these directions, U,(x) is the external flow, and r = r(x) 
denotes the distance of a point on the boundary from the axis of symmetry, 
as in Fig. 1. 

The equations proposed in (1) which transform the coordinates and 
velocities of the axially symmetrical problem to those of the equivalent 
two-dimensional problem are 


(5) 


(6) 


(7) 

U, = U,, (8) 
where L denotes a constant length. 

Consider now the problem of a wall jet on a body of revolution, formed 
when a circular jet impinges upon it along its axis. Equation (8) is satisfied 
since U, = U, = 0 as there is no external flow and hence the two-dimen- 
sional analogue, independent of r(x), is the plane wall jet investigated 
in (2). The solution for the plane wall jet, given in (2), is 


(9) 


(0) 


| r(x) da, | 
0 
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Fic. 2. The velocity function f’(y). 
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with the stream function % given by 


= (7). (11) 

The arbitrary reference velocity U may be expressed in terms of the flux 

of exterior momentum flux, which is invariant for any given two-dimen- 
sional wall jet. Thus if F is given by 


al at ag] dj, (12) 
é 
then U = 40F/v*. The function f(#) satisfies the equation 


f" +I" +2f = 0, (13) 
and the boundary conditions 


f'(~)=9. (14) 
It is shown in (2) that the solution with f(oo) = 1 may be selected without 
loss of generality. The velocity function f’(y) is shown in Fig. 2. 

We may now apply the transformation given by equations (5)-(8) to 
calculate the characteristics of a wall jet on any body of revolution, as 
the one shown, for example, in Fig. 1. Equations (6) and (12) show that 
for such a wall jet the flux of exterior momentum flux is still constant, 
and thus we may define a constant F which is related to this quantity by 


F= | ay| dy. (15) 


Equation (15) may be deduced directly from equations (3) and (4) by 
multiplying (3) by r and integrating from y to 00, then multiplying the 
resulting equation by ru and integrating between 0 and oo. 

Before proceeding further it is convenient to modify equation (5) thus: 
| r(x) dx +C, (16) 
0 
where C is a constant length. The constant C may be interpreted as a 
change of origin in the two-dimensional flow, which is obviously allowable 
physically; but since there is no external flow the introduction of C makes 
it possible to incorporate an analogous arbitrary length in the solution 
for the axially symmetrical case. 

Taking now as a particular example the plane radial wall-jet, with 
r(x) = x, the results of (2) and (3) may be recovered. If the constant 
length L is given by 


t= 75 


v 
(17) 


j 
| | 
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we have, from (6) and (16), 1 ae a constant length, 


= (18) 


9 = (19) 
where C = U*/3/3v?, and so from (7), (9), and (10) 


¥ 


27vat 
The arbitrary reference velocity U may be determined, from (15), in terms 
of the flux of exterior momentum flux as 


(22) 


_ 18F 
=| (23) 


135F2* 

which are exactly the results obtained in (3) by very different methods 
(the definition of F in (3) differs by a factor p® from that given here, p being 
the density). The length / is of the order of the distance taken by the 
wall jet to attain its final similarity form. Incorporating the length / in the 
solution means that the velocity given by (23) remains finite as x > 0 and 
so may have some sign.ficance even near the axis; a method whereby / may 
be roughly determined is given in (3). 

As a further example we may take as the body of revolution a hemi- 
spherical shell of radius a with r(x) = asinx/a. The solution applies to 
either a wall jet formed on the outside or the inside of the shell, assuming 
that the jet adheres to the surface; for if the flow separates from the surface 
the results no longer hold. We then have 


[2x —asin(22/a)+ (25) 


j= ysin(xja), (26) 


10F 
10F sin‘(x/a) ] 


(27) 


(28) 


201 
C = — 
40F’ 
giving 
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where / is a constant length as before. Equations (27) and (28) indicate 
the variation of velocity and wall jet thickness with z. 

The characteristic properties of wall jets on any body of revolution can 
be similarly calculated. It is interesting to note that, apart from a change 
of scale, the velocity profile is the same for all such bodies of revolution. 
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A THEORETICAL STUDY OF THE DIFFUSION OF 
TRACER GAS IN AN AIRWAY 
By D. W. JORDAN 


(Mining Research Establishment, National Coal Board, 
Worton Hall, Isleworth) 


[Received 21 July 1960] 


SUMMARY 

The theory of G. I. Taylor relating to diffusion in a turbulent stream is used to 
provide constants appropriate to the diffusion of gas in an airway. The diffusion 
equation is solved for steady and instantaneous gas sources of various types in 
airways of arbitrary cross-section, and the mixing length for a source waved about 
in the airway is estimated. The work relates to the measurement of air flow in 
airways when the use of anemometers is not possible. 


1. Introduction 


A METHOD in use for estimating the flow of air in mine roadways where the 
use of anemometers is inappropriate is to release at a point in the airway 
a steady or instantaneous supply of a suitable tracer gas (1). This mixes 
with the air stream by turbulent diffusion, and at a distance downstream 
where mixing is complete, measurement of concentration in the case of 
steady emission, or of the integral of concentration with respect to time in 
the case of instantaneous emission, provides a measure of the mean velo- 
city, or of the volume flow of air. 

If the emitter provides tracer at a steady rate e per unit time and the 
concentration is measured at any point at a sufficiently large distance 
downstream, the mean air velocity @ is given by 


= (1.1) 
where c is the concentration and s the cross-sectional area. If an instan- 


taneous source emits an amount m, the concentration-time curve c(t) is 
noted and @ is given by 


m = sé { c(t) dt. (1.2) 


The distance required to ensure uniform mixing, however, is often incon- 
veniently large, and this study evaluates various possibilities for reducing 
this length. Among the methods considered are the best position of the 
emitter (or receiver), the use of an extended source (or receiver), and the 
production of good initial mixing by waving the source about in the air 
stream. 


(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 2, 1961] 
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The theory of turbulent diffusion in pipes is given by Taylor (2). Taylor 
takes the cross-sectional distribution of stream velocity and turbulent 
diffusion coefficient appropriate to turbulent pipe flow, shows that the 
dispersion of concentration in the axial direction is largely governed by the 
variation of stream velocity over the cross-section, which produces a 
greater spread axially than would be expected from turbulent diffusion 
alone. The net effect of turbulent diffusion and axial convection is, at a 
sufficiently large distance, approximately the same as if the concentration 
were being dispersed axially by a diffusion process only, described by a 
coefficient about 100 times as large as the coefficient of turbulent diffusion 
[see also Aris (3)]. 

In the present work we assume generally that the diffusion of tracer gas 
in the air is described by a constant radial coefficient k,, and a constant 
axial coefficient k,, which is much greater than k,. Numerical values are 
assigned to these constants at a later stage. It is realized that it is doubtful 
whether the Taylor theory applies accurately at the moderate distances 
downstream which are considered, but there is evidence (4) that the 
accuracy is sufficient to support the trends noted in this paper. 

Firstly, some rosults are proved for the general case, where the stream 
velocity and diffusion coefficients are assumed to vary on the cross-section, 
and later a group of particular problems are solved on the assumption that 
mean values of these quantities may be used as an approximation. 


2. General results, with non-uniform stream velocity and diffusiv- 

ity and arbitrary cross-section 

2.1. Flux, emission, and concentration 

Let the concentration at the point (x, y,z) be c(z, y,z,t) per unit volume, 
the stream velocity in the z-direction be v(x, y), the coefficient of turbulent 
diffusion in the (x, y)-plane be k,(x, y), and in the z-direction k,(x,y). The 
equation governing the variation of concentration is easily shown to be 


a) Af, a) af, a) 


in regions free from sources, and the boundary condition implying no loss 
at the wall is 


ec 
k 0, 21.2 
( ) 


where the derivative is normal to the wall. 

Two points of some importance will now be considered. Assume firstly 
that a steady emission is produced by any distribution of sources. Then 
by integrating (2.1.1) over the cross-section (noting that the time deriva- 
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tive is zero), and applying Green’s theorem and the boundary condition 


(2.1.2), we have dady = 0, 


in any region which does not contain sources, where the integral is over 
the area of the cross-section. The integral is equal to the flux f of diffusing 
material over the section. We have then 


0. 


(2.1.3) 


Therefore downstream of all sources the flux has a constant value f+ and 
upstream another value f~-, and 

ft-f-=e, (2.1.4) 
where ¢ is the total emission of the sources. Methods of measuring velocity 
assume that f- is zero, though for very small air speeds some loss to —90 
might be expected. But for this region 


f-=[[ re. dedy, 
where c_,, is the concentration as z-> —oo. But since v is positive and 
there are no sources at z = —0oo, c_,, and therefore f~, must be zero. The 


downstream flux can therefore always be equated with the total emission 
in the steady case. 

To estimate the flux, equations (1.1) and (1.2) are used, and their 
accuracy depends upon measurements being made far enough downstream 
for the concentration to have become uniform. We note that the errors 
involved in the two methods due to relaxation of this condition are pre- 
cisely the same. Defining m,,, the (erroneous) estimation of the integrated 
mass crossing the section 


= 00 f edt, 
0 


where the instantaneous source is emitted at ¢ = 0 and the measurement 
is made at the fixed point (2, y,z), then integrating equation (2.1.1) with 
respect to time and observing that 


0, 


we see that m,., obeys the steady-state equation of diffusion. The boundary 
condition on m,,, is also identical with the steady case. The two methods 
are therefore precisely equivalent under all circumstances, and are also 
equivalent to a method involving an emission lasting a finite time. No 
further reference will be made, therefore, to the method described by 
equation (1.2). 


> 
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2.2. Distribution of sources required to produce a uniform concentration 
at all points downstream 
We consider the function of space and time defining a ‘unit instanta- 

neous point source’ emitted at time ¢ = ¢, at the point (29, ¥o,2)). For 

t < t, the function is zero everywhere, and for t > t, it satisfies equation 

(2.1.1) and the boundary conditions. It is continuous with continuous 

derivatives, except that, for t = t), at the point (79, Yo, t)) it has a singu- 

larity expressing the fact that the total quantity of diffusing substance 
released into the medium is unity. We now prove the following theorem: 

A unit instantaneous point source at the point A on the cross-section z = Z 
released at t = ty, produces at the point B on the cross-section z = z, at time 

t = t,, the same concentration as if the positions of A and B on their respec- 

tive cross-sections were interchanged. 

This result is already familiar in the heat conduction case, where the 
velocity is zero and the diffusion coefficient uniform (5). If 


G (a4, Yr, ty | Lor Yor Zo» te) 


represents the concentration at y,,2,,t,) due to a source at (29, Yo, Zo, 9) 
the result to be proved can be written 


G(X, Yas by | Los Yor 20 to) = Yor — 29; | —%1, (2.2.1) 
Writing (5) 
y,z,t) = G(x, y,z,t | Xo, Yq, to), zero for t < ty 
G(x, y, z,t) = G(x, y, —2z, —t|2,y,, —2z,, —t,), zero for t > 
the equation satisfied by G, is 
eG, 2 eG, 
= (2.2.3) 
where 8(s) is the Dirac delta function. G, satisfies the adjoint equation 


== (2.2.4) 
It is easily seen that 


G, L(G)—G, Z(G,) 


which proves the statement regarding adjointness. This in turn is equal 


5: 

= 
= 

. 
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to the expression obtained by applying the same process to the right-hand 
side of (2.2.3) and (2.2.4): 
(2.2.6) 
By integrating (2.2.5) and (2.2.6) over all space and time, and using 
Green’s theorem to transform the second and third terms of (2.2.5) to 


a surface integral (which has value zero by the boundary condition 
(2.1.2) satisfied by G, and G,), we obtain 


= 21 — Yo» Zo (2.2.7) 
The left-hand side of this equation is zero, since G, is zero for t = —oo 
and G, is zero for t = +00, and both are zero for z = +00. The right- 
hand side is therefore also zero, and by (2.2.2), this is the result (2.2.1) 
which was to be proved. The slight integration difficulty regarding the 
delta functions can be avoided by regarding them as the limit of one of 
their analytical approximation functions with continuous derivatives; for 


example, 8(8) as A> oo. 


The result can be extended to the case of steady emission from a point 
source or, by integration, to sources and receivers which are distributed over 
a finite area or volume. Instantaneous sources and continuous receivers 
cannot, however, be interchanged. 

The symmetry of the Green’s function, proved above, can be used to 
calculate what cross-sectional distribution of steady source strength is 
needed to produce a uniform cloud at all points downstream. 

The concentration arising from a distribution of sources of strength 
€(Xp, Yo) (emission per unit time), on the plane z = 2, is given by 


e(x,y,2) = (xo, | Yor%) (2.2.8) 


where y, Z | %g, Yq; Zp) is the steadily-emitting counterpart of the instan- é 
taneous point source above. The value of the function ¢€(79,¥9) which % 
makes c independent of (x,y,z) is suggested by the relation defining the a 
flux: 
ft= If dady = constant (2.2.9) 


If we identify «(x», with the operator 


¥ 
ay 
i 

a 

>. 
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where C is a constant, and apply the reciprocal relation to equation 
(2.2.8), we have 


e(z,y,z) = C | Lor Yor Zp) 
0 


0 


where Zo=—% and Z= —z. (2.2.12) 


v(Zo, Yo)—ky Yo x, Z) dx, dyy, (2.2.1 1) 
0 


But by equation (2.2.9), this is proportional to the flux over the plane 
Z, from a unit steady source at (x,y, Z), which is unity, multiplied by C. 
The function ¢(x, y,z) is therefore equal to C, and is independent of (x, y, z). 

The first and second terms in equation (2.2.10) refer respectively to a 
distribution of sources of strength proportional to the stream velocity at 
the point, and to a distribution of dipoles of strength proportional to the 
axial diffusion coefficient at the point. The dipole term could generally 
be ignored owing to its smallness at fairly short distances downstream. 
A similar result applies for instantaneous source distributions, the quantity 


y, 2, t) dt, (2.2.13) 
0 


integrated from the instant of emission, being independent of x, y, and z. 


3. A duct of arbitrary cross-section, with physical parameters 
supposed independent of position on the cross-section 
Analytical solutions of the general equations of the last section for initial 

value problems are not easily found. Numerical solutions are possible, but 

difficult to apply to variants such as a moving source. In the following 
sections it is assumed that no important qualitative errors will arise by 
treating the velocity and diffusion coefficients as having constant values 
over the cross-section, the axial diffusion coefficient being modified to take 
into account differential axial convective transfer by accepting as a value 

Taylor’s ‘virtual diffusion coefficient’. The merits of these approximations, 

inasmuch as they concern the radial distribution at moderate distances, 

are something of an open question: in the case of a steady source in a 

circular duct, reasonably good agreement of the theory with experiments 

(made at the Mining Research Establishment of the National Coal Board) 

has been found, but for instantaneous sources there is no evidence available. 
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In a circular duct the mean radial diffusion coefficient (2) is given by 

k, = 0-04avyt, (3.1) 
where a is the radius, v the mean stream velocity, and y a dimensionless 
coefficient depending on the roughness of the walls (2). Taylor's ‘virtual 


diffusion coefficient’ for combined diffusion and differential convection in 
the axial direction, for fully developed turbulent flow, is 


ky = (3.2) 


3.1. Point sources in a duct of arbitrary cross-section 


If x and y are rectangular coordinates in the cross-section, and z is 
measured in the direction of flow, the equation of diffusion can be written 
dc 


where 


It is convenient to define dimensionless parameters as follows: 
T = (k,/k,)*vt/l, Z = (k,/k,)t2/l, X = Y =y/l 
K = (k,k,)t/vl, C = Pe, (3.1.2) 


where / is some typical length associated with the problem, for example, 
the radius in the case of a circular airway. Equation (3.1.1) then becomes 


ec 12aC 
p= 0 (3.1.3) 


where V? now means é7/0X?+-é?/aY, and 
W = Z-T, 


(3.1.4) 


and the boundary condition is 


—— = 0, (3.1.5) 
where @/én denotes the derivative normal to the bounding surface. 

The equation for an instantaneous point source of unit emission in these 
dimensionless variables (emission (k,/k,)* in the dimensional variables), at 
the point W = 0,X = X,Y = Y, is 


1 


P 


a q 
“aa 
3 
Be: 
50092 .54 
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Applying a combined Laplace and Fourier cosine transformation, and 
defining 


(X,Y, £,p) | coséW dW | e-O(X,Y,W,T)dT, (3.1.7) 


the resulting equation is 
p/K}O = (3.1.8) 


The solution to this equation is the Green’s function of the Helmholtz 
equation, which can be expanded in terms of the eigenfunctions of the 
homogeneous Helmholtz equation as follows (5). The homogeneous 
equation 1s V2p-+-Ay = 0, (3.1.9) 
with the boundary condition 


(3.1.10) 


This equation has solutions only for a real, discrete set of values of A; 
A = 0 is clearly an eigenvalue, with the corresponding eigenfunction equal 
to A, a constant. The others are numbered n = 1, 2,..., 00, in order of 
increasing magnitude. The corresponding eigenfunctions can be arranged 
to form an orthonormal set, such that, integrating over the cross-section, 
lifm=n 
= 11 

Assume that this has been done. Then this requirement gives 
A= S+, (3.1.12) 


where S is the area. The expansion of the Green’s function in terms of 
the eigenfunctions is then 


pe 2 l (X » Yo, An Win(X, ¥,A,) 


where in the present problem, 
—(é2+p/K). (3.1.14) 
C can now be recovered from C by applying the transformation which 
is the inverse of (3.1.7). Equation (3.1.13) is in a particularly suitable 
form for the application of the calculus of residues. The terms in (3.1.13) 
have singularities at 4? — A2, giving simple poles at 
p = —K(€*+A§). (3.1.15) 
We apply the inverse transform for the parameter p and evaluate the 
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residues from (3.1.15) in the usual way. The inversion of the resulting 
expression with respect to é can be got from tables of cosine transforms (6). 
The final result for the concentration arising from an instantaneous point 
source is 


TYAKT 
~ “Ont ist > 


The concentration arising from a point source emitting continuously at 

a unit rate is easily obtained by integration with respect to the time, or 

by omitting the time variable in the above analysis. The resulting 
expression is 


(This corresponds to a rate of emission equal to v/?(k,/k,)* in the original 
dimensional system.) 


An approximate expression, equivalent to (3.1.17), is obtainable at 
moderate distances downstream. If 


Z/K = vz/k,> 1, (3.1.18) 


and the distance downstream is such that only the first few terms of 
(3.1.17) are required, so that in addition 


4.2 K? = 42k, /lv <1, 


eb Zl 
Yn(Xo, Yo, (3.1.17) 


(3.1.17) becomes 


and this is also the expression obtained by neglecting in the differential 
equation the term @C/éW?, i.e. neglecting axial diffusion in the steady 
case. The exponent in (3.1.19) is independent of k,. This indicates that 
in the steady case the concentration should not be sensitive to the value 
assumed for the axial diffusion coefficient. 

For calculations in particular cases the eigenfunctions and eigenvalues 
for use in the above equations can be obtained by analytical or numerical 
methods. At moderate distances downstream, which is the case of practical 
interest, only the first or second will be required because of the rapid con- 
vergence of the series. 

It should be noted that, in the case of a steady source [equation (3.1.17) ] 
the mean value of the concentration over all cross-sections is equal to its 
value at Z = o. 


n=1 q 
4 
F 
| 
a 
‘ 
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3.2. A duct of rectangular cross-section 
We require the eigenfunctions and eigenvalues of equation (3.1.9), for 
the rectangular region 


0<X<A, O0<Y<B, (3.2.1) 
satisfying (3.1.10). The normalized eigenfunctions are 
$,, = 2(4 B)-* / B), (3.2.2) 
PP 


where p and q are positive integers, not both zero. The expression for an 
instantaneous point source in the plane Z = 0 at the point (Xp, Y,), 
emitted at time 7' = 0, is 

~ SA * 
(3.2.4) 
and the approximate expression for a steady point source (equation 
(3.1.19)) is 


x 14+4> > cos( prX pr X/A)cos(qx¥ /B)}. 
(3.2.5) 


In the case of a square cross-section each eigenfunction is counted twice, 
since 


A2(p,q) = AX(q, p). 


3.3. Point sources in a circular duct 
We require normalized eigenfunctions of the Helmholtz equation with 
the boundary condition (3.1.10) at R = 1 which are finite in the cross- 
section and are also periodic in the polar angle. These functions are 
To(%om 2) Jol 
R, 8, Ay) = | 2tcosné J, R)/7*( — 
24 sin NO (3-3-1) 


Here we have A, = ad, (3.3.2) 


and «,,, is defined by 
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The concentration from an instantaneous point source at (Ry, 4), on 
the cross-section Z = 0, is given by 


eZ-TYAKT 


J, (Xam R)J, R ) 


and the approximate expression for a continuously emitting unit source is 


n=1 >O 


Aom>O 


For a point source placed centrally, Ry = 0, and the double sum in 
(3.3.5) vanishes, giving 


n=1 


Finally, an expression for a ring source concentric with the circular 
duct is obtained. 
A continuous distribution of instantaneous point sources along R = Rg, 
0 < @ < 2m is equivalent to an instantaneous ring source, and the result 
is easily obtained by integrating equation (3.3.4), and is seen to be the 
same as if all the terms for » > 1 were simply omitted. Similarly the 
approximate expression for a continuously emitting ring source of total 
emission unity (in the dimensionless units) is, from (3.3.5), 


Xom>0 
When R, = 0 this reduces to the result for a central point source. The 
reciprocal theorem proved above for the general case is obviously applied 
in all these expressions. 

A check on the theory for a steady, centrally placed point source was 
provided by experiments performed at the Mining Research Establish- 
ment of the National Coal Board (S. Shuttleworth, unpublished work), 
in a smooth-walled duct, with nitrous oxide as tracer gas. The duct was 
one foot in diameter, the airflow 13-37 ft/sec, and the estimated value 
for the resistance coefficient y, considering the duct to be smooth, was 
0-005. k, can be obtained from equation (3.1). The results are shown in 


7 
ae 
| 
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Fig. 1, and it can be seen that reasonable agreement was obtained with 
the theory. The approximate expression (3.3.6) was used. The fact that 
the experimental points generally lie somewhat closer to the mean is 
perhaps due to the fact that the duct was not perfectly smooth. 
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Fic. 1. Comparison of theory with experiment for a steady central point source in 
a duct of diameter | ft, showing concentration profiles at various distances down- 
stream. Experimental points denoted by rings. 


When a central point source is used the distance required to obtain 
mixing adequate for the use of equations (1.1) and (1.2) is often greater 
than is convenient under mining conditions, The use of concentric ring 
sources can improve this situation considerably. In Fig. 2 are shown the 
concentration profiles arising from ring sources of dimensionless radii R,, 
obtained from equation (3.3.7), at various distances downstream. In the 
case R, = 0 (equivalent to a central point source, the curves all cross at 
a point about 0-63 radii from the centre, where the concentration is equal 
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to the mean concentration. This has been observed experimentally, both 
in the experiments described above, and in work by Hodkinson (7). It 
arises from the fact that the first eigenfunction J,(a), R) of equation 
(3.3.6) vanishes when R = 0-627, so that only the subordinate terms in 


(bt) R= 0-4 


AR, =O 
(CENTRAL POINT SOURCE) ad 


o-2s (4) 


on T T T T 
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Fic. 2. Theoretical concentration profiles, expressed as ratios (concentration/mean 
concentration), of ring sources radius Ry. The numbers on the curves are values 
of K’Z’ = 0-04yiz/a. 


the series give a contribution to the deviation from the mean concentra- 
tion. By the reciprocal theorem above, a ring source of radius equal to 
this number should give effectively uniform concentration at all points 
downstream, and this tendency is exhibited in Fig. 2 (c), where Ry = 0-6. 
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4. Acceleration of mixing by using a moving source 

The above theory shows that the concentration arising from a steady 
point source does not become uniform and equal to its mean value over 
the cross-section until a fairly large distance downstream. By moving the 
source about in the airway in a random manner it might be expected that 
a good degree of mixing would be obtained near the plane of the source 
and the final mixing thereby accelerated. The question is investigated in 
detail for a circular airway, but the results apply qualitatively to an 
airway of any cross-section. 


4.1. Concentration due to a moving source 

Suppose a source in the plane Z = 0 is moved about whilst emitting at 
a constant unit rate, its position at time 7' = J, being given in polar 
= Ry = R(T) (4.1.1) 
and the process started at 7’ = —oo. The concentration due to an equiva- 
lent source over the interval d7,, at the point (#, 0, Z, 7’) is 

C{R, Z, T—T, Ro(T)} 

where C is the function (3.3.4); and the concentration arising from the 
moving source is 


T T 
Cn = CAT | (4.1.2) 
where 
€, exp[ —(Z—T)?/4KT —a2,, KT] 


and l ifn=0 


2 ifn>0 
(4.1.3) 
The term in (4.1.2) under the summation sign represents the fluctuating 
part of the concentration, which it is desired to minimize. Consideration 
will, for simplicity, be confined to the prospect of improving the evenness 
of the cloud in a region sufficiently far downstream so that only the first 
term in the summation is of importance. In the case where the receiver 
can be located at any point, the lowest eigenvalues are «,, and a,,, but 
since the receiver would normally be expected to be placed at or near the 
duct centre, the terms which involve these eigenvalues vanish at this point 
and the next larger value, giving the principal term, is 
a = a, = 3-832. (4.1.4) 


3 
— 
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At the centre we may therefore write the fluctuating part as 


T 
(1) = aR) KT] 


f(T) = JofaRo(T)} 


4.2. Scanning the cross-section by the source 

To ensure that the average value of the function W(7') should be zero, 
it is easily proved to be sufficient that the mean value of f(7,) should be 
zero. If a path R,(7)), 0(7,) is selected, which covers all elements of cross- 
sectional area to an equal extent, then the mean value of f (T,) will be zero. 
This is true, not only for the single term of (4.1.2) selected above, but for 
all the terms, and indeed for the corresponding terms in a duct of any 
cross-section, for 


y, Ay) dady = ds = 0, 


by Green’s theorem, using the boundary condition on the eigenfunctions. 
The details of the fluctuations will be confined to the single case above, 
though the general argument will apply in all cases. 

It is convenient, for simplicity, to consider that the scanning by the 
source is performed by means of a series of random cycles of equal period, 
each of which attempts, however imperfectly, to cover the whole area of 
cross-section. The cycles could be chosen to cover the area densely, in 
which case the mean value f(7,) over a cycle would be very nearly zero, 
though each cycle would take a long time: or very sketchy cycles could 
be adopted, each being performed quickly, so that /(7j) over a cycle 
would not be very close to zero, but the mean over many cycles would 
approach zero. A typical cycle of each type is shown in Fig. 3. Cycles 
of the latter type only [Fig. 3 (6)] will be considered here. 

It can be proved that, at the moderate distances which we are consider- 
ing, the function g(7) is very nearly a Gaussian curve [Fig. 3 (c)] with 


standard deviation op, = (2KZ)\. 
Its ‘width’ may therefore be measured by 
AT, = 2e7, = 2(2KZ)}. (4.2.1) 


(This is equal to the ‘duration’ of a cloud from an instantaneous point 
source.) It can be seen from the integral (4.1.5), that if a fairly large 
number of random cycles, each of duration 57,, are completed in the 
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4. Acceleration of mixing by using a moving source 

The above theory shows that the concentration arising from a steady 
point source does not become uniform and equal to its mean value over 
the cross-section until a fairly large distance downstream. By moving the 
source about in the airway in a random manner it might be expected that 
a good degree of mixing would be obtained near the plane of the source 
and the final mixing thereby accelerated. The question is investigated in 
detail for a circular airway, but the results apply qualitatively to an 
airway of any cross-section. 


4.1. Concentration due to a moving source 

Suppose a source in the plane Z = 0 is moved about whilst emitting at 
a constant unit rate, its position at time 7T' = 7, being given in polar 
and the process started at 7’ = —oo. The concentration due to an equiva- 
lent source over the interval d7,, at the point (R, 6, Z, 7’) is 

C{R, Z, Ro(T)} 

where C is the function (3.3.4); and the concentration arising from the 
moving source is 


T T 
Cm = CAT =-+ 47, (4.1.2) 


where 
Snub To) R,(T,)}cos 
and _ {1 ifn=0 
\2 ifn>9 


(4.1.3) 


The term in (4.1.2) under the summation sign represents the fluctuating 
part of the concentration, which it is desired to minimize. Consideration 
will, for simplicity, be confined to the prospect of improving the evenness 
of the cloud in a region sufficiently far downstream so that only the first 
term in the summation is of importance. In the case where the receiver 
can be located at any point, the lowest eigenvalues are a,, and a,,, but 
since the receiver would normally be expected to be placed at or near the 
duct centre, the terms which involve these eigenvalues vanish at this point 
and the next larger value, giving the principal term, is 


= = 3832. (4.1.4) 
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At the centre we may therefore write the fluctuating part as 


T 
where (7) — aR) T— 
out Thx) (KT) (4.1.6) 


= T,)} 


4.2. Scanning the cross-section by the source 

To ensure that the average value of the function W(7) should be zero, 
it is easily proved to be sufficient that the mean value of f (7) should be 
zero. If a path Ry(T)), 0(7Z,) is selected, which covers all elements of cross- 
sectional area to an equal extent, then the mean value of f (7,) will be zero. 
This is true, not only for the single term of (4.1.2) selected above, but for 
all the terms, and indeed for the corresponding terms in a duct of any 
cross-section, for 


by Green’s theorem, using the boundary condition on the eigenfunctions. 
The details of the fluctuations will be confined to the single case above, 
though the general argument will apply in all cases. 

It is convenient, for simplicity, to consider that the scanning by the 
source is performed by means of a series of random cycles of equal period, 
each of which attempts, however imperfectly, to cover the whole area of 
cross-section. The cycles could be chosen to cover the area densely, in 
which case the mean value {(7}) over a cycle would be very nearly zero, 
though each cycle would take a long time: or very sketchy cycles could 
be adopted, each being performed quickly, so that /(7) over a cycle 
would not be very close to zero, but the mean over many cycles would 
approach zero. A typical cycle of each type is shown in Fig. 3. Cycles 
of the latter type only [Fig. 3 (6)] will be considered here. 

It can be proved that, at the moderate distances which we are consider- 
ing, the function g(7)) is very nearly a Gaussian curve [Fig. 3 (c)] with 


standard deviation on, = (2KZ)}. 
Its ‘width’ may therefore be measured by 
AT, = 207, = 2(2KZ)t. (4.2.1) 


(This is equal to the ‘duration’ of a cloud from an instantaneous point 
source.) It can be seen from the integral (4.1.5), that if a fairly large 
number of random cycles, each of duration 57), are completed in the 
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duration, AT, of the function g(7’—7,), and if the value of f(%) over 
each cycle is not very different from zero, then the value of W(T7') will 
also be very nearly zero, independently of the time. This is a more 


Track of source: 


(a) 


Tracks of source in each cycle 


Fic. 3. (a) Values of f(7,) during a long cycle with dense area coverage; (b) values 
of f(T) over a succession of short cycles of period 57,; (c) the shape of the curve 
g(T —T,), of width AT. 


favourable condition than that its mean value should be zero. On the 
assumption that 8T, < AT, (4.2.2) 


we can write approximately, using 57, as the integration interval, 
T 


W(T) = | AT = (4.2.3) 


T, 
T, 
: 
(c) 
9(T-T,) 
r T, 
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where the summation is carried out over all the cycles, [9(77—7)] is the 
mean value of g over the sth interval 57), and 


(a) 
taken over the corresponding cycle duration 57). It can be shown that 
provided the scanning speed is constant 


| dT, = (4.2.4) 

(3) 
where H, is a quantity depending only on the path of the source during 
the sth eycle, and not on its duration (see Appendix). Now suppose that 
the values H,, taken from a population consisting of a great many random 
and independent cycles (this was the original hypothesis for the motion 
of the source), have the mean value zero (which will be the case if the 
track is random with respect to the area elements); and that they are 
distributed, not necessarily normally, about this value with variance oj. 
Then random values of the function W(7') will have mean value zero, 
and a variance, ojy, given by 


wv => 57, |. (4.2.5) 
This sum may be recast as an integral if (4.2.2) holds, so that 


T 
The integral can be evaluated explicitly, giving at R = 0 (the centre) 
4.2.7 


The spread of values of W(7') is then measured by 2c,,. We wish to 
measure the improvement on the case when the source is stationary; say, 
at the centre of the duct. This measure is the ratio of 20, to the first 
term under the summation in equation (3.3.6), which gives the deviation 
from the mean due to the stationary source. We finally have: 


2x standard deviation of concentration 
‘improvement _ from mean value (moving source) 
ratio’ " steady deviation at the sampling 
position (fixed central source) 


= 2(8T,)to,, *K-*Z-* = (4.2.8) 
0. 
by equation (4.2.1). 
In order that the above analysis should be applicable, 57,/A7, must be 
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small. It seems likely that a value < } should be sufficiently small to 
ensure that the arguments hold, and in this case 
improvement ratio < 0-6loy, (4.2.9) 
and the smaller this quantity the greater the improvement over a fixed 
source. The requirement (4.2.2) makes it difficult to ensure that adequate 
cycles, giving small values of o,, can be achieved in the time available, 57). 
It is therefore necessary to consider the actual time limit imposed by 
(4.2.2). The expression for AJ), (equation (4.2.1)) in ordinary time units, is 
At, = 2-8Atatyt (4.2.10) 
where y is the Chezy—Darcy resistance coefficient, mentioned above, and 
A is a constant which for fully developed turbulent flow has Taylor’s 
value 7-14. When the flow is not fully developed, A is somewhat smaller 
(4), and to take a very unfavourable case, we put 
A = 2-0. 
Also, under mining conditions, it is known that y lies generally between 
0-06 for rough, obstructed airways, and 0-01 for smooth airways, so that 
equation (4.2.10) becomes 
At, > 1-95(z/a)'v/a_ for rough airways 
> 1:25(z/a)tv/a for smooth airways 
The requirement that 57,/AT7, should be less than 4 implies 
57, < 0-33(z/a)tv/a_ for rough airways 
< 0-21(z/a)tv/a for smooth airways 
Thus, for sampling at 20 diameters downstream in a 10-ft diameter 


airway, with a velocity of 3 ft/sec, the time available for a cycle is about 
2 sec at most. In Fig. 4 are shown some typical cycles, each consisting 


H= -O9 ‘03 ‘10 
Fic. 4. Typical cycles of motion of the source, with corresponding values of H, 
the measure of the coverage achieved. 


of only three straight line movements of the emitter, which could easily 
be accomplished in this time, with the associated values of H, (see Ap- 
pendix). These values are small, of the order of 0-1, and it is to be expected 
that, for such cycles, oy = O(0-1). (4.2.11) 
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Using this value, equation (4.2.8) gives 
improvement ratio = O{0-061} = O{6%}, (4.2.12) 
so that with even such poor coverage as is provided by these movements 
of the emitter, the spread of concentration at the centre of the duct about 


its mean value is only of order 6 of that due to a steady source at the 
centre of the duct. 
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APPENDIX 
THE CALCULATION OF H, 
Over a typical cycle 


To+8To To+ OT 
Te To 
Suppose the emitter is waved at constant speed V, and that the total length of 
its track in a cycle is L, = V8T, (A2) 


then (Al) may be written 
Le 
0 


where L(T,) is the length covered between T' = 0 and 7%, and is an increasing 
function. This may be rewritten 


1 
R=0 
where dz, L is the total length of the path between R, and R,+85R,. If the whole 
area were covered evenly with a dense sequence of tracks, which is the ideal case, 
then 


dp, L = 2L,R,dR,, (A4) 
making (A3) zero. In the case of imperfect coverage we have 
dp, L = R,+ x(R,)} (A5) 
where x measures the deviation from the ideal case. Putting then 
L,V- = 87, (A6) 
Tot+8To 1 


we have dT, = AR, = 87, H,, (A7) 
0 


Af 
‘ 
5 
4 
4 
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where H, is a quantity measuring the adequacy with which the track of the sth 
cycle covers the area, and is zero for a perfect track. It is independent of time. 
To estimate the values of H, for the tracks in Fig. 4, we put 


$( Ry) = (A8) 


where L( R,) is the total length of path inside the radius R,, and, integrating (A7) by 
parts: 


i 
H, = (Re) + Ala). (A9) 
0 


¢(R,) is a function that can easily be measured from diagrams. The reader is 
reminded that a is the first zero of the equation 


Liz, 
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ON THE APPROXIMATE CALCULATION OF HEAT 
TRANSFER BY LAMINAR FLOW FROM A ROTATING 
HEATED DISK AT PRANDTL NUMBERS BETWEEN 
0-2 AND 6 


By D. R. DAVIES (University College of North Wales, Bangor) 
and C. B. BAXTER (University of Sheffield) 


[Received 19 November 1959. Revise received 24 October 1960] 


SUMMARY 

Calculations of heat transfer from rotating disks and surfaces of revolution in 
practical circumstances are likely to involve surface temperature distributions which 
are functions of radial and meridional distance from the axis of rotation. Exact 
solutions are not available for these boundary conditions, and the available approxi- 
mate methods have been shown (1) to break down at values of the Prandtl number 
o less than 6. In this paper a new approximate method is presented which is shown 
to be accurate (maximum error 5 per cent) for values of o between 0-3 and 3 (this 
range includes the important case of air, for which o = 0-72) and a reasonably 
accurate one (maximum error 10 per cent) for values of o between 3 and 6. It is 
applicable for a fairly general form of radial surface temperature distribution, leading 
to a convenient expression for the local rate of heat transfer, and can probably be 


extended to include the calculation of heat transfer from rotating surfaces of 
revolution. 


1. Iv has been shown by Davies (1) that a linear, tangent form of approxi- 
mation to the exact radial velocity component on a rotating disk leads 
to accurate predictions of the rate of heat transfer for values of o (the 
Prandt) number) greater than 200, and a power law approximation leads 
to reasonably accurate predictions for values of o between 6 and 200. 
These results are due to the fact that for very large values of o the ratio 
of the thickness of the temperature boundary layer to that of the velocity 
boundary layer is small; the crucial part of the latter is very near the 
surface and is sufficiently well represented by the linear approximation. 
For smaller values of o it is necessary to obtain a good representation of 
a greater proportion of the velocity profile and this is obtained by a power 
law approximation. However, for values of o less than 6 (including the 
important case of air, for which o = 0-72), the part of the velocity profile 
which lies farthest from the surface of the disk becomes of increasing 
importance. In this region the radial velocity, having reached a maximum 
value from zero at the surface, decreases to zero again in the surrounding 
stationary fluid. A power law representation can no longer be used but 
we find that, by multiplying the power law by a suitable expression, a 
[Quart, Journ. Mech. and Applied Math., Vol. XIV, Pt. 2, 1961) 
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profile can be constructed which rises to a maximum and decreases again 
to zero, and is a good representation of the exact profile. An analytic 
solution for the temperature distribution, leading to a convenient expres- 
sion for the local rate of heat transfer, can then be obtained for the case 
when the radial surface temperature distribution can be represented by 


the fairly general form 3 a,, 7”, where r denotes distances from the centre 
=0 
of the disk. js 


2. If z denotes distance measured normally from the surface of the 
disk, u, the velocity component in the radial direction, u, the velocity 
component normal to the surface, the stream function, v the kinematic 
viscosity, x the thermal diffusivity, 7 the temperature, the temperature 
equation in von Mises form is 


eT @ 

—_ = — 1 

(1) 
where the frictional heating term is neglected, and the boundary layer 


approximation applied. 
An exact calculation of u,, by Cochran (2), is available and is expressed 


u,/(wr) = F(2), (2) 


where { = (w/v)*z; w is the angular velocity of the disk. The problems 
involved in treating equations (1) and (2) by electronic computation 
methods for general boundary conditions are not yet solved and, in order 
to obtain solutions, we require an approximate representation of the exact 
values of u,/(wr) in terms of %; this must be accurate over the central core 
of the velocity profile for values of o of the order of 1-0, and reasonably 
accurate over the whole width of the profile. 
We therefore extend previous work by Davies (1) and write 


u,/(wr) = C ( | (3) 


where the parameters B, C, t, and s are selected so that (i) an analytic 
solution of the temperature equation is obtainable, (ii) the continuity 
relation is satisfied, and (iii) a good fit to the exact profile is obtained. 

In order to satisfy the first of these requirements we substitute (3) into 
(1), which becomes 


(4) 


where == fils, t+s=2, and D = 


CALCULATION OF HEAT TRANSFER BY LAMINAR FLOW 225 


The velocity approximation, represented by (3), must satisfy the con- 
tinuity relation in the form 


ru, = (5) 
and substituting (3) into (5), and integrating, leads to 


f dp = CBN, 


where = = (6) 


The form (3) approximates closely to a power law at small values of C. 
This suggests we select values of C and t close to those used previously by 
Davies (1); we find such suitable values are t = 0-4 and C = 0-4. An 
appropriate value of B can quickly be calculated by first differentiating 


(3) in the form u,/(wr) = CB(1—¢2-). (7) 


This shows that u,/(wr) reaches a maximum at ¢ = (2/t)"?-, Substituting 
this value of ¢ and the exact maximum value of u,/(wr) into (7) yields the 
value B = 0-64. The integral in (6) is then evaluated numerically at 
suitable intervals in ¢; it is convergent as ¢ > 0, since t < 1, but as ¢— 1, 
it increases indefinitely so that as ¢ > 1 (u, > 0), { > 0, as in the exact 
calculation by Cochran (2). Finally, by expressing the values of ¢ in (7) 
in terms of {, we obtain an approximate velocity profile. 


TABLE 1 


Exact values of u,/(wr), calculated by Cochran, and approximate values, 
calculated from equation (7) with t = 0-4, B = 0-64, and C = 0-4 


© 025 O§ OFF 10 125 11°95 BO 235 25 30 39:25 


Exact © O100 0-153 0178 O172 0137 O116 O-100 0°084 0°074 0061 0-050 
Approx. | © O-180 O175 O59 O°137 0°097 0°079 0°067 0°052 


This profile is expressed numerically in Table | and is seen to be a good 
approximation over almost the whole range of ¢, the maximum deviation 
from the exact profile being not more than 4 per cent for values of ¢ 
between about 0-25 and 2-25. This is clearly the crucial range of {, when 
o is of order one, and the temperature field extends over the whole velocity 
field. 


3. Analytic solutions of (4) are most conveniently derived by separation 

of variables for the fairly general case in which T,—7, = s a,,1”, T, being 
n=0 

the temperature of the surrounding stationary fluid and 7, that of the 


5002 .54 Q 


0 . 
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disk surface. We write 7’'—T, = 0, so that the boundary conditions become 
=0, when §=1, and @= a,7, when (8) 
n=0 


By separating the variables in (4), we see that a suitable solution of 
the form 


a, (a, 1/8; €)+E, 2—1/s; £)] 
a= 

(9) 
exists, where F is a hypergeometric function [see, for example, Sneddon 
(3)|; the parameters « and f are given by a+8 = 1/s—2Ds and a8 = nD, 
and HE, = —F(a,B; 1/8; 1)/F(a+1—1/s,B+1—1/s; 2—1/s; 1), the hyper- 
geometric series in the numerator and denominator of this expression 
being convergent (3) for all values of o and n and rapidly convergent for 
the value of o for air (0-72). The local rate of heat transfer Q at distance 
r from the centre is then given simply by 


(s—1)CB'-*k(w/v)t E,a,r". 


Conversely, if the rate of local heat transfer Q is prescribed and can be 
represented conveniently by a power series, the coefficients a, can be 
determined to give the associated temperature distribution. 

It is, however, difficult to construct solutions by this method when @ is 
any other function of r on € = 0, but formal solutions can be obtained by 
means of Laplace transforms, when @ is a function of the form g(r)—g(0) 
along € = 0, g(r) being any continuous function. 


4. We suppose, for mathematical convenience, that heating of the disk 

commences at r = ry, say, and that along £ = 0, we have 
= T—T, = g(r)—9(r0) 

for r > ry, and @ = 0 for r < ry; this restriction to the form g(r)—g(r,) 
ensures the removal of a discontinuity at r = rp. 

We then write log(r/r,) = R, so that R varies from 0 to «0, and if we 
multiply each term in (4) by e~?”, and integrate with respect to R from 
0 to co, we obtain the ordinary differential equation 


a6 fi dé 
where §=— | de-PR dR. 
0 


The transformed boundary conditions are 
§6=0 on §=1, 


= 
4 
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and §= | dR = f(p,r,) on &=0. (11) 
A suitable solution is 


= f(p,r9)[ F(a, B; 1/8; F(«+1—1/s, B+-1—1/s; 2—1/s; €)], 


where now af = pD; applying the inversion theorem [see, for example, 
Carslaw and Jaeger (4)], we then obtain the formal solution 


ytio 
| BE F(a+-1—1/8, B+-1—1/s; 2-1/8; £)] x 
In the case when @ = $ a,r along ¢ = 0, this solution can be shown to 
n=0 


approach the form given by equation (9), as ry > 0. 


5. An exact solution is only available in the particular case when the 
surface temperature is independent of r, and in order to obtain an estimate 
of the accuracy of our approximate method, we compare the approximate 
and the exact numerical values in this case. We write 


where 7; is the constant temperature on the surface of the disk, so that 
6é=1o0n &=0 and @=0 on €=1. A solution of (4) exists which is 
independent of r: this is 


é 
= f ag 2Ds), (13) 
0 


where A(1—1/s,2Ds) is a beta function [see, for example, Sneddon (3)]. 
The total rate of heat transfer Z for a disk of radius r is then given by 


= 2(2—t)C[(1—t) B'“B(1—1/s, 2Ds)]-2. (14) 


This expression has been computed for various values of o using t = 0-4, 
B = 0-64, C = 0-4, and the results are shown in Table 2. The correspond- 
ing exact values for EZ have been given by Millsaps and Pohlhausen (5) 
for values of o between } and 10; they are easily calculated (when frictional 
heating is neglected) in this case for values of o less than 4, and values 
down to o = 0-1 are also given in Table 2. 

Comparison of the results shown in the table indicates that our approxi- 
mate method is accurate (the maximum error is 5 per cent at o = 3) for 
values of o between 0-3 and 3 (this range includes the important case of 
air, for which o = 0-72) and a reasonably accurate one (maximum error 
is 10 per cent at o = 6) for values of o between 3 and 6. 


(7 —T)/(T,—T) = 8, 
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TABLE 2 


Calculated values of rate of heat transfer from a rotating disk 
using the exact and approximate solutions 


Calculated values of E/{k(T,—T,)eotv—4r2] in cm. sec units 
Approximate values (using Percentage underestimate 
Prandtl (14) with t = o-4, B = 0-64, of the exact values by the 
number Exact values C = 0-4) approximate formula 
or 0°35 0-30 14 
0°54 8 
0°67 3 
o's 0°86 0°86 ° 
I 1-26 1°24 I 
2 1°66 1-62 2 
3 2-04 1-92 5 
4 2°36 2°16 8 
5 2°62 2°36 9 
6 2°84 2°55 10 


Finally, we note that our method of approximation can probably be 
extended to calculate heat transfer by laminar (and possibly turbulent) 
flow from rotating surfaces of revolution for Prandtl numbers of order one, 
when almost the whole velocity profile has to be taken into account and 
the tangent (Lighthill) and power law (Davies) methods of approximation 
break down. We note also that the effect of frictional heating can be 
added to the present solution [see, for example, Millsaps and Pohlhausen 
(5)]. 
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THE STEADY-STATE RESPONSE OF CONTINUOUS 
DISSIPATIVE SYSTEMS WITH SMALL 
NON-LINEARITY 


PART I 


By R. D. MILNE (Queen Mary College, University of London) 
[Received 24 May 1960. Revise received 29 September 1960] 


SUMMARY 


A theoretical method is developed for calculating the periodic response of con- 


tinuous, dissipative systems which exhibit a small degree of non-linearity in the 
stress-strain relation. 


1. Introduction 


In recent years it has become increasingly important in the field of 
aeroplane dynamics to be able to estimate both the nature and distribution 
of damping forces throughout metal structures. In particular, the two 
important problems of flutter (oscillatory instability involving aero- 
dynamic forces) and noise-induced vibrations require such information 
not solely for the introduction of dissipation as such but also for the 
change in the mode of motion that may ensue. 

The success of linear vibration analyses of structures, both simple and 
built-up, together with the fact that the damping is always small in 
magnitude, leads to the conclusion that a theory based on a small devia- 
tion from the generalized Hooke’s law will probably suffice to describe 
the behaviour of metal structures. Thus in this paper the stress at a point 
is taken to be the sum of a linear and a non-linear function of the strain 
at the point, the non-linear contribution being small. 

Various forms of stress-strain law have been proposed to describe the 
dissipative forces in metals: a comprehensive account of the probable 
mechanisms of these is given in (1). For a general (transient) motion 
these laws generally include terms representing hereditary effects. How- 
ever, if the motion is restricted to be steady (in the oscillatory sense), 
these hereditary effects are also steady and the cyclic stress is effectively 
an ‘instantaneous’ function of the cyclic strain. The analysis of this paper 
is restricted to steady motion. 

Early tests (2) on structural metals indicated that the damping forces 
were proportional to the amplitude of motion but in phase with the 
velocity. These results have formed the basis for a linear theory of 
(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 2, 1961) 
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structural damping usually called Hysteretic Damping (3, 4, 5) since the 
stress~strain locus describes a closed curve (an ellipse) during each cycle 
of the motion. The area enclosed by the curve is proportional to the energy 
dissipated per cycle. 

In this paper the form of the non-linear stress as a function of strain 
is not restricted in any way (except in the assumption of steady motion) 
but the main interest is fixed on stress-strain curves which may be double- 
valued. In the event that the form of this curve is independent of the 
frequency it may be determined statically. Any general double-valued 
curve may lead to non-linear stiffness in addition to non-linear damping 
forces. 

Experiments on built-up structures (6, 7) show a marked non-linearity 
in damping and this is usually attributed to the presence of riveted or 
bolted joints: idealized models of such joints may be dealt with by the 
analysis of this paper. 

The presence of non-linearity means that the response to a harmonic 
force is periodic but not harmonic. This analysis is based on an expansion 
of the response as a Fourier series the coefficients of which are functions 
of position. These are determined using the so-called averaging principle of 
Ritz, a method which has been used by Klotter and others (8, 9) for the 
analysis of dynamical systems having non-linear elements. 

It may perhaps be emphasized that in the analysis of continuous bodies 
the mode of response is of vital interest, particularly in those problems 
in which aerodynamic pressures on the body themselves depend on the 
motion of the surface of the body. 

In Part II the analysis is applied to the extension, torsion, and flexure 
of slender bars, and a numerical example is presented for the case of 
torsion of a bar with riveted joints. 


2. Analysis 
(i) Series solution for periodic response 

Referred to rectangular cartesian axes O(2,,2,,2,) and in the absence 
of body forces, the equations of motion for an element of a continuous 
body are 


00 
pti, —# — 0 (j,k = 1,2,3), (2.1) 
Ox; 


where p is the density of the body, u,; the displacement vector, and o;, 
the stress tensor and the summation convention is adopted. The displace- 
ment is small so that the strain tensor is given by the linear relations 


Ou, 
= +—}. 2.2 


: 
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For a virtual variation 5u; in the path of the motion the total work 
done by the forces acting on every element vanishes at every instant; 
integration over any arbitrary time interval t, < t < t, yields the varia- 
tional equation of motion 


V 


in which the volume integral is taken throughout the body. 
An application of the divergence theorem transforms (2.3) to 
ts 
ptt, du; dV + on dey, dV— bu ds} dt=0, (2.4) 
i, 
wherein s; = o;,l,, the 1, being the direction cosines of the normal to the 
body surface. 

It will often be convenient in the sequel to represent the six independent 
stress and strain components by ¢,, e, respectively: the summation con- 
vention will apply with r, s = 1,..., 6, and 

= 1,2,3)=ey = &), 
(r = 4,5,6) = 2ey, (j #k). (2.5) 

Throughout the analysis it is assumed that a sufficient number of points 
of the body are fixed so that a rigid-body motion is not possible: this 
assumption is made on the grounds of convenience only. 

Let the surface tractions s; have simple harmonic time dependence of 
circular frequency w but arbitrary phase and magnitude over the surface. 


Using a bar superscript to denote the amplitude of a simple harmonic 
quantity, then 


8; = §,coswt = (5+ 15/)cos wt. (2.6) 
The linear operator ¢ is defined by the relations 
isin wt = cos wt, icos wt = —sin wt, (2.7) 
so that in particular, xi=#= 


and the operator has the characteristic values -+,/(—1): consequently the 
quantity (a+b) may be treated algebraically as a complex number. 

The important restriction is imposed that the solid under consideration 
is ‘essentially dissipative’. That is, in the absence of external forces, any 
motion, whatever its amplitude or form, will eventually decay; the possi- 
bility of internal self-excitation is thus expressly excluded. With this 
restriction a steady-state response to harmonic external forces will be 
established after a sufficient time which will itself have a definite periodicity. 
Every variable then describes a closed path in the phase plane and in 


4 
— 

k 


232 R. D. MILNE 


particular time integrals representing hereditary effects do not appear 
in the description of the motion. The periodic strain being known at 
every point throughout the cycle it follows that the stress may be uniquely 
specified at every point of the cycle; thus while the stress may be a double- 
valued function of the strain as discussed in section | it is a single-valued 
function of time within the fundamental period. 

The stress-strain relation will be represented by 

o, Crs (2.8) 
where the C,, are constants and the functions f, are known either for all e, 
(single-valued) or only when further conditions on the strains are specified 
(double-valued) such as e, positive or negative, increasing or decreasing in 
magnitude. The deviation from the generalized Hooke’s law is measured 
by the parameter » which will be considered to be so small that yu? and 
higher powers of 1 may be neglected compared with unity. Obviously the 
actual value of » is dependent on the absolute values of the strains in 
any particular case: it is here assumed that C,,e, and f,(e,) are of the same 
order of magnitude. 

The analysis proceeds on the basis that the steady-state response to 
external excitation of circular frequency w will have the fundamental 
period 7’ = 22/w. Subharmonics of the forcing frequency are at present 
excluded: this restriction is discussed later. 

Accordingly, the response u,(t,t+-7') is represented in the interval 
t, <t<+t,+T by the series 

= 008 aut = (iat itl; at. (2.9) 

The independent displacement vectors @;,, a; are determined by 
satisfying the variational equation (2.4) approximately by choosing the 
hitherto arbitrary time interval (¢,—t,) to be the periodic time 7’. Since 
(2.4) must be satisfied for every virtual variation 57a, p We obtain the set 
of (time independent) variational equations 


— (Pw)? J pil, dV+ | dV — | = 0, 
V 
(2.10)t 


in which the surface integral is non-zero for 8 = 1 only. The complex 
stress F_, is given by 


2 é,, 008 awt }cos Bust d(wt) — 
0 


2a 


+ An asterisk denotes the complex conjugate. 
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The stress Fg may be expressed in terms of the strains @,, by the use 
of complex operators defined for steady periodic motion as represented 
by (2.9). Thus (2.10) can be written in the form 


— pit, dV + J (Onn +Draglng — J [5 dS = 0. 


(2.12) 
The complex operators D,,2 = K,,+iG,,g are given by an integration by 
parts of (2.11); thus 


2a 
(2.13 a) 


-- 
Al 008 08 Bt d(wt). (2.13b) 
The individual coefficients K,,.9, G,,g are clearly not uniquely defined by 
the use of the derivative form, but since the equations of motion in- 
variably involve only the summation D,,,é, this particular synthesis is 
nov important in itself. The importance of the forms (2.13) lies in the 
fact that they render the actual calculation of the non-linear contribution 
from the stress functions f, much easier than if (2.11) is used. 


In (2.13) write cos ast = 1,,_008 a(wt-+¢,.) (2.14) 
and change the variable to 
The expression for K,,9, for example, then becomes 


where f,, denotes éf,/ée,. Substituting 


2a 
= J + 2 008 x 


2 008 a), g}sin df, 2, (2.164) 


and similarly, 


= 2 A,qp 008 x 


x 


= 

2a 
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A, = COS —itsin 
= 1, , +i sin ad, ,)(cos ag, sin ag, g) 
,(cos Bd, sin Bd, 

The D,,g are thus functions only of the A, ,g and it will be noted in 
particular that when a = then A, 9, = The case a = = 1 leads 
to the results obtained (for systems with a single degree of freedom) 
by application of the method of equivalent linearization of Kryloff and 
Bogoliuboff (10). 

The practical interpretation of the formulae (2.16) in some particular 
cases is discussed and illustrated in section 3 of Part IT. 

The Dg are, of course, functions of position, implicitly through the 
spatial variation of the r,, and ¢, , and explicitly through any variation 
there may be in the defining functions f, throughout the body. 

The variational equations (2.12) imply the following series of differential 
equations and boundary conditions for the di, g: 


(Bw)* pit; ba, = 


(B= 1) 917 
The usual boundary conditions will specify a; , = 0 over a region S’ of S 
for all 8 with &,; specified over the remaining region (S—S’) of S. 


(ii) Form of the response in particular cases 

The general problem described above is clearly extremely formidable 
and an obvious simplification would be to terminate the series (2.9) after 
an arbitrary, small number of terms. However, by considering the probable 
physical behaviour of the system by analogy with its linear basic system 
(« = 0) and utilizing the condition of small non-linearity, a number of 
more firmly based approximations emerge which lead to considerable 
simplification. 

Consider the basic linear system belonging to the body under con- 
sideration and described by the differential equation of motion 


+-— (Cir fim) = 9, 
ex, 


with a, = 0 on S = S’ and zero surface traction. This equation describes 
the simple harmonic motion of the linear body of circular frequency A 
under no external forces: it possesses an enumerable infinite set of non- 
trivial solutions a, ,, (the eigenfunctions) for particular discrete values of 
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the parameter A* = A? (the eigenvalues). The eigenfunctions are the normal 
modes of the body. Under external harmonic excitation of frequency Q, 
the motion (of frequency Q only) is completely specified by a weighted 
sum of the eigenfunctions, where the weighting numbers are given by 


Q, 
— 
in which M, = J dV, = ! dS. 


As 02 > 8, q, > 0, and a; > a@,,. Further, consider the response of the 
special linear system described by the equation 


+ +4) Fim} = 9, 
where » is a small constant. In this special case the eigenfunction expan- 
sion still holds but with 


and in particular, when Q? = 2, 


= 
Vn M,, 
In this case it has been assumed that the linear dissipative force is distri- 
buted in exactly the same manner as the elastic stiffness: when this is 


not the case the simple eigenfunction expansion fails, but near the un- 
damped resonances we may still assert that 


O(1/p). 


x O(1/p). 


In this connexion the deductions made by Rayleigh (11) concerning the 
motion of linear dissipative systems having a finite number of degrees of 
freedom are relevant. 

The following arguments concerning the non-linear system are based 
on the foregoing results in conjunction with the assertion that the system 
contains dissipative forces of O(,). 

We distinguish the following cases for the response of the non-linear 
system described by (2.17) to external harmonic forces of frequency w. 


(a) Non-resonant response 

None of the eigenvalues of the linear basic system approach the 
forcing frequency w» or any of its multiples, that is, A,, ~ Bw. In this case 
the fundamental response @,;, to O(1) surface forces is O(1): possible 


: 
3 
O(force) 
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‘excitation’ forces in any harmonic d,, represented by the non-linear 
coupling terms D,,,é,, are thus, at most, O(u) and the response a, ¢ 
is O(u). Furthermore, since resonance in the fundamental is excluded, 
an expansion of @,, in powers of is valid so that for the total response 
we may write (to O(,)), . 


u,(t,t+T) = 008 008 (2.18) 
where U, is O(1). Similarly, writing 
e,(t,t+7') = E,,008 wt +p 008 (2.19) 
then by a Taylor expansion (2.8) becomes 
0, = B, 008 > é,,..008 +pf,(B, cos wt). 


Since the expansion in powers of y is valid for all y the set of differential 
equations (2.17) becomes 


purl; +5 E,,,} =0 (2.20 a) 
with unit order surface tractions §; and 
+Fing( Es, = 0 (8 = 1,2,3,...) (2.20b) 


with zero surface traction. 

The unit order component U; is thus the solution of the linear basic 
system and the O(.) displacements @,;, are determined from the (linear) 
equations (2.20b) which are fully specified whenever 0; is known. This 
type of power expansion was introduced by Lindstedt (12) for the solu- 
tion of non-linear systems in free motion: in this context it is invalid 
whenever natural frequencies of the linear system are approached by the 
forcing frequency. 


(b) Resonance in fundamental only 

Let w be close to an eigenvalue A,,; then the response in the fundamental 
to external force of O(u) will be O(1). Provided Bw + A2, (8 = 2,3, 4...) 
then the harmonic response @,;, will be the non-resonant response to 
excitation forces of possible O(j) and will therefore be O(u). Thus in this 
case the harmonics are all O(u) times the fundamental and the appropriate 
forms of (2.9) and (2.8) are 


u; = %;, coswt+p > a; , cos awt, (2.21) 
0, = 008 wt +p aut} coset). (2.22) 


| 
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The differential equations (2.17) become 


s+ + Pima} = © (2.23 a) 
with an surface traction and 


(B = 2,3,4,...) (2.23 b) 


with zero surface traction. The equations (2.23b) are fully specified on 
solving (2.23 a). 


(c) Resonance in fundamental and several harmonics 

Let w be close to an eigenvalue A,, but such that f’w = A,, also for one 
or more values of Bf’, m. Then the response to O() external forces is O(1) 
in the fundamental and also in those harmonics for which 8 = f’; the 
appropriate form of (2.9) is then 


U; = 2 COS a’ wt +p 2 COS (2.24) 


The differential equations (2.17) divide into the two sets 


é 
p(B + + Preamp = 9 (2.252) 
with O() surface forces for 8’ = 1 only and 


+ + 8) Fim =0 (B48) (2.25b) 
with zero surface forces. 


(d) Resonance in harmonics only 

In this case w 4A, but p’w =A,, for one or more values of f’, m. 
Fundamental response to O(1) forces and harmonic response for 8 = f’ 
are O(1) so that this case is in effect identical with case (c) in so far as (2.25) 
hold but the order of magnitude of the response is lower. 

It may be remarked that, in the foregoing, whenever dissipative forces 
are dependent on the time rate of change of strain, the forces D,, 2é, ¢ 
contain the quantities Bw as factors, with the result that the harmonics are 
effectively more heavily damped than the fundamental: case (6) may then 
often be a sufficient approximation to problems strictly falling within 
case (Cc). 

The possibility of resonance occurring simultaneously in the funda- 
mental and harmonics suggests the possibility of response occurring near 
the resonances of the linear system lying below the excitation frequency. 


| 
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In this case a result similar to that of sub-harmonic resonance (12) in 
systems with a single degree of freedom would occur in which the response 
to excitation of period 7 has periodicity y7', where y is an integer. Thus, if 
an eigenvalue of the linear system approximates to w/y the expansion 
(2.9) for u; may be modified to 


with the same analysis as before. The presence of more than one frequency 
coincidence, A,, = w/y, leads to rapid lengthening of the gravest period, 
since it must be commensurate with all the ‘sub-harmonics’ assumed 
present. 

There is no fundamental objection to the further development of the 
analysis to include excitation by non-harmonic periodic forces which may 
be represented by Fourier expansions: a number of the equations (2.12) 
then contain a forcing term. 

The foregoing arguments may be extended to this case but with less 
simplicity. Non-linear interaction between harmonics will still be governed 
in a general way by the rules already applied. In particular, whenever 
the linear system has resonant frequencies which are integral multiples, 
then strong coupling between harmonics may be expected. 


3. Self-excitation by external forces 

The analysis thus far has been developed in terms of the steady-state 
response of the system to an external harmonic force whose magnitude 
is independent of the motion of the system. The presence of external 
forces which are linearly dependent on the (surface) motion will not alter 
the analysis in principle, provided the motion is steady. In the absence 
of any independent external force the system is one in ‘free’ motion and 
the restriction to steady motion limits the discussion to those conditions 
of the system for which a steady motion is possible, that is, a condition of 
neutral (oscillatory) stability. Sufficient disposable parameters must be 
available to satisfy this condition. Alteration of these from their ‘critical’ 
values will lead to oscillatory stability or instability. 

The extremely important aeroelastic phenomenon of fiutter falls into 
this category; here the aerodynamic forces are dependent on the surface 
motion of, say, a wing oscillating in a uniform stream of air. The aero- 
dynamic theory is generally linearized, and furthermore the unsteady 
aerodynamic forces are usually calculated for simple harmonic motion 
only (hereditary effects due to the downwash field of the wing are steady 
in the oscillatory sense). In the linear case the stability boundaries of the 


| 
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motion are invariably examined by assuming the motion to be simple 
harmonic and then determining the simultaneous values of the frequency 
and air speed which render this motion possible (double eigenvalue prob- 
lem) together with the mode shape. An increase of speed above this 
value will lead to instability. This method can readily be applied to the 
case of small non-linearity discussed in this paper. The system acts as 
a selective filter so that only the flutter mode and frequency are significant 
in the feed-back loop from body motion to surface forces. Thus for the 
investigation of flutter all terms in the expansion (2.9) may be omitted 
except one which we may conveniently take to be the fundamental; that 
is, u,(t,t+T7') = d;, coswt. Equations (2.12) and (2.17) apply with 8 = 1 
only and with &; wholly dependent on the surface motion. 

It may be noted that, in fact, the aerodynamic forces are valid only for 
the single term solution. 

Flutter of dynamical systems with mechanical units having small non- 
linearity is dealt with in (13). 


4. Solution by the Rayleigh-Ritz procedure 

A very common approximate method of solution of problems derived 
from a variational principle and the one used almost universally in aero- 
elastic work is the Rayleigh—Ritz procedure in which the spatial displace- 
ment field is constrained to be represented by a weighted sum of chosen 
coordinate functions. 

Since each differential equation (2.17) has the same linear part and 
satisfies the same kinematic boundary conditions there seems no reason 
for not expressing each harmonic displacement function a; in terms of 
the same set of spatial coordinate functions v;,,. Thus, we set 


Pa 


in which the complex numbers a,, represent the phase and modulus of 
the real coordinate functions v;,,. The suffix p, on the summation sign is 
to indicate that the summation may contain only some of the set v, , for 
any given harmonic. Since the a, , are independent and da, , 4 0 each 
of the variational equations (2.12) leads to a set of algebraic equations, 
the typical fth set having the form 


is J = 0, (4.2) 


‘9 : 
2 
f 
j 
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where the e,, are the strain components corresponding to the displace- 
ment fields v,,. In terms of generalized coefficients the fth set of equations 


where 
Myq = J AV, dogg = J esp 
= aq Q = J 2S. (4.3) 


It will often be convenient to use a number of normal modes of the 
linear system as the coordinate functions v,,, and in this case m,, = k,, = 0 
for p and k,, = 

The D,,g are of course functions of the r,,, ¢,, (equation (2.13)) and 
through the relations 


ap, +t sin abs.) > (4.4) 
Pp 


are implicit functions of position and of the generalized coordinates a,, ,. 

The relation (4.4) is clearly not sufficient to determine the r, , and ¢, , 
in terms of the a,,, except when the right-hand side contains one term 
only. This means that an iterative solution is essential beginning with a 
single term approximation in each harmonic: this single term approxima- 
tion must then be used to determine the D,, , for subsequent approximations 
containing more terms in the coordinate function series. An example 
illustrating this procedure is given in section 4 of Part II. 

The approach in the flutter problem will be somewhat different since 
(a) only the fundamental is considered, and (6) the frequency is to be 
determined. Taking as many coordinate functions », ,, in the fundamental 
as seems advisable the equations may be solved for the linear case » = 0, 
and the frequency, air speed, and flutter-mode shape found. Using this 
mode and specifying the absolute magnitude of one generalized coordinate, 
say @,,, the D,, may be calculated. The solution is now found to this 
‘linear’ problem. Keeping the magnitude of a,, fixed at the previously 
chosen value the flutter mode may be improved upon by successive solu- 
tions. When the mode has stabilized the specification of a,, implies that 
the absolute magnitude of the displacement at every point is known. 
Repeating the complete calculation for different values of a,, a plot of 
flutter speed and frequency against displacement amplitude will be the 
final result. 

Both the forced and flutter solutions outlined above are suitable for 
computation on a digital computer, the successive approximations being 
obtained by the solution of sets of linear equations. 


|| 
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THE STEADY-STATE RESPONSE OF CONTINUOUS 
DISSIPATIVE SYSTEMS WITH SMALL NON-LINEARITY 


PART II 


By R. D. MILNE (Queen Mary College, University of London) 
[Received 24 May 1960. Revise received 29 September 1960] 


SUMMARY 


The analysis of Part I is applied to the motion of slender bars and a numerical 
example is discussed in some detail. 

In addition the practical application of the general method of analysis is briefly 
discussed. 


1. Introduction 


Ir is shown in Part I (1) that if the displacement vector at a point is 
represented by the series 


u(t,t+7) = > &, , cos awt (1.1) 


(1, equation (2.9)), then the harmonic displacement vectors @,, satisfy 
a coupled set of non-linear differential equations. The non-linearity is 
embodied in the complex operators D,,, whereby the stress components 
in the Sth harmonic are given by the expressions (C,,+D,,2)é,, where 
the C,, are functions of position only (the linear elastic constants) and 
the D,,2 are functions of position and the strain harmonics é,,. By the 
assumption of small non-linearity < C,,. 
The D,, are calculated from the non-linear part of the assumed stress— 
strain relation through the equations (2.13) of Part I. 


2. Application of the analysis to slender bars 
The general analysis of Part I is adapted to the small motion of slender, 
isotropic bars for which the total displacement field is the sum of suitable dis- 
placement fields representing longitudinal extension, bending, and torsion. 
Consider a uniform bar having as coordinate axis Oz, the line of centroids 
of the cross-sections, with the ends of the bar at z = 0,1: the axes Ox, Oy 
are taken to be parallel to the principal axes of the cross-section. 
The following system of displacements is adequate to describe the 
deformation of the bar when the rotation of any element is small (2): 
u = U(z)—y6(z) 
v = V(z)+20(z) 


(Quart. Journ. Mech. and Applied Math., Vol. XIV, Pt. 2, 1961] 


(2.1) 


ais 
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where @ is the angle of rotation of sections about Oz and ¢ is the Saint- 
Venant torsion function for the cross-section (this function ensures com- 
patibility only for a linear stress-strain relation so that its use here is 
approximate). We may set V(z) = 0 for simplicity. 

It is convenient to introduce the flexural axis passing through the point 


1 1 
"jj 
where P, = yd, P, dA 
A 


by writing (2.1) as 
Up—(y—yr)?, 
v= 


dV, 
The torsion function ¢, = ¢—zy,+ yz, satisfies the relations 


= tes = drydA = 0. 
A A 


The displacement field (2.1) leads to the three non-zero strains 
_ dW d*Up d*6 


= 2e,, = + (2.2) 
= = Ox —(y 


For periodic motion the expansion (1.1) is applied to the displacement 
systems of (2.1) so that, for example, 


6= > cosawt, ete. (2.3) 
a 


Substituting in the variational equation (2.12) of (1) for the Bth har- 
monic with = 80$—(y—yr) 56%, etc., 


we obtain three differential equations and their associated boundary 
conditions for the three displacement systems Ug, Ws, and 4g. 

In deriving these equations the simplifying assumption is made that 
the density p is uniform over the cross-section: in this way a large number 
of (in this context) unnecessary inertial coupling terms are omitted. The 


4 

‘ 

3 

| 

| 
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surface work integral (1, equation (2.12)) takes the form 


i 
0 


where X, Z, © are the components of external loading and the external 
torque per unit length respectively: this external loading is non-zero for 
8 = lonly. The quantities Ng, 8, Ms, Tz are respectively the end load, 
shear component, bending moment, and torque on the bar: since these 
may depend on end conditions involving the total motion of the bar they 
may appear in each harmonic. 

The fth set of differential equations and associated boundary conditions 
are: 


(a) In Ug, 


with 
dU, dU, d dW, 
aU, dW; d | 
(b) In Wa, 
(2.5) 
aW, dW, 
(c) In 8g, 


dz?| dz Pdzt' dz 
d aw; aU, dé 
+ P88 + | =0 (2.6) 
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where @ is the angle of rotation of sections about Oz and ¢ is the Saint- 
Venant torsion function for the cross-section (this function ensures com- 
patibility only for a linear stress-strain relation so that its use here is 
approximate). We may set V(z) = 0 for simplicity. 

It is convenient to introduce the flexural axis passing through the point 


| waa, Ur = | 


where P, = dA, 
A 


by writing (2.1) as 
Up—(y—yr)?, 
v = 


wl tur dV, 
The torsion function ¢, = ¢—zy,+ yx, satisfies the relations 
= [ dA = dry dA = 0. 
A A A 
The displacement field (2.1) leads to the three non-zero strains 


For periodic motion the expansion (1.1) is applied to the displacement 
systems of (2.1) so that, for example, 


Up = cos awt, 
6= > 6, cosawt, ete. (2.3) 


Substituting in the variational equation (2.12) of (1) for the Sth har- 
monic with bus 86%, 


we obtain three differential equations and their associated boundary 
conditions for the three displacement systems Ug, Wg, and 9g. 

In deriving these equations the simplifying assumption is made that 
the density p is uniform over the cross-section: in this way a large number 
of (in this context) unnecessary inertial coupling terms are omitted. The 


C3 = Cx, 
| . 2.2 


= 


STEADY-STATE RESPONSE OF DISSIPATIVE SYSTEMS. II 245 
surface work integral (1, equation (2.12)) takes the form 


where X, Z, © are the components of external loading and the external 
torque per unit length respectively: this external loading is non-zero for 
8 = lonly. The quantities Ng, 8.2, M,,, Tz are respectively the end load, 
shear component, bending moment, and torque on the bar: since these 
may depend on end conditions involving the total motion of the bar they 
may appear in each harmonic. 

The fth set of differential equations and associated boundary conditions 
are: 


(a) In Us, 
dW, aU, 
with 
dU, @U, d 
p, dW, d 
(b) In Wp, 
(2.5) 
dW, dW, 
with [#4 “e+ (Pep d = 
(c) In 8p, 


dW, aU, a6 dé 


dW, 


PU, dé 
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with 
dO, d 


By 
| — Bator as 


In the above equations ZH, G denote the Young and Shear Modulus 
respectively, A is the area of the cross-section, and [and J are the torsion- 
bending constant and torsional rigidity given by 


| $444, 
A 
J= J le + + | dA. 


The symbols ®,, , are all functions of the displacement systems Us, 8g, Wa, 
and are defined by the following integrals taken over the cross-section: 
= [ Dug dA (n = 0,1, 2), 
4 


O55 


me l,n= 0 


= | a" dA with =Ia=1 
A 


Op m = 2,n = 0, 


A 
m=0,n=0 
with 
m= i, 9 = 0, 


Das| dA, (2.7) 


The equations (2.2) show that any harmonic strain component @, , is 
the sum of products of a real function of z, y, say f(x,y), and a function, 
say g(z), of z which is complex in the phase vector sense; thus for a typical 


t 
= 
If \g,| is the modulus of the vector g then, clearly, 


ap, +i sin ap 4) = f(x, Y)|9q\(cos ag,-+é sin (2.8) 
and thus ry = f(x,y) \gq(2)|- 
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Thus in the D,,g we may replace the r, in the A,g by the products 
f(x, y)\g,| and understand the ¢,(z) to refer to the phase of the vector 
g(z). The ®, , then take the form 


= F(x, f(x,y) |Ga(2)|,ba(2)} 4A, (2.9) 


and carrying out this integration gives 


These general results are illustrated in the example of section 4. 

The differential equations show that, even in the absence of inertial 
coupling (2» = yy = 0), bending and torsional motions may be coupled 
through the non-linear terms: the coupling with longitudinal motion is 
practically of less importance. Coupling arises through: 

(a) possible dependence of direct stress on shear strain and vice versa 
(see, for example, quadratic and cubic stress-strain relations in classical 
elasticity (2), 

(b) through direct stress terms involving warping. 

When end constraint against axial displacement is not important, then 
the term in e, (equation (2.2)), which is proportional to rate of twist, may 
be omitted and then only 9,3, ®, 3, and are non-zero. 

It is common to apply the displacement field equation (2.1) to non- 
uniform beams provided the variation of cross-section is nowhere abrupt: 
the equations require little modification to meet this case. 

It will often be the case that the non-linear stress functions are them- 
selves functions of position in the cross-section in addition to being 
dependent on local strain amplitude; that is, the non-linear ‘cross-section’ 
of the bar will differ from the normal linear ‘cross-section’, Thus, in 
bending, the ratio of the non-linear to the linear forces in the cross-section 
is represented by ®,/HP,, while in torsion the corresponding ratios are 
®,/GJ and ®,/ET. It follows, for example, that the dissipative forces in 
bending vibrations arise mainly from the material situated farthest from 
the neutral axis of bending. 

Equations (2.4), (2.5), (2.6) are extremely involved even when we 
assume uncoupled bending and torsional motions; hence, it is necessary 
to resort to a Rayleigh—Ritz procedure except in the simplest cases, As 
an example, consider the form of the equations for the case of uncoupled 
bending vibrations when we also neglect the small effect of rotatory inertia 
(the term pP,»(d*U,/dz*) in (2.4)). That is, we set 


= 2 (2.10) 


t 
: 
: 
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where the v,(z) satisfy at least the kinematic boundary conditions at 
z = 0, 1; the equations of motion are then given by (1 (4.3)) with the 
modified definitions 


= | pAv, v, dz, diag = dz, 
0 


i I 
= BP. = | Xo, ds (B=1). (2.11) 
0 


3. The calculation of the non-linear stress coefficients 
(a) Fourier analysis of the stress 

The results obtained by the method of section 2 (1) based on approxi- 
mately satisfying the variational equation of motion may be reproduced 
by a direct expansion of the non-linear stress f, as a Fourier series. Thus, if 


uf, = F,,..008 awt, (3.1) 


then the Fourier coefficients F,g are given in the usual way by the formula 


Qn 
= { f,008 Pet det) —i J st 
0 


which is in accord with (2.11) of (1). 

The stress, however, is directly dependent on the strain and it is thus 
natural to refer a harmonic expansion of the stress to the phase and 
amplitude of the corresponding harmonics of the strain: the coefficients 
of this series then play the role of operators which so multiply and rotate 
the strain vectors that the appropriate stress components are obtained. 
To this end set 

(3.2) 
so that (3.1) becomes 


Differentiate both sides of this equation with respect to wt and equate 
the sth partial derivative on the left side to the sth term on the right, thus 


of, de, 


d 
dt) = 2, Feat) 008 cunt) 


= — ,(K,,,, 8in awt+G,, , cos at). 
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Multiplying both sides by é, sin Bwt and integrating over the interval (0, 27), 
we obtain 


ef, de, i =— 


which is in accord with formula (2.16a) of (1). Similarly, multiplication 
by é, cos Bwt and integration yields the formula (2.16 b) of (1). 

The way in which the stress harmonics are related to the phase and 
amplitude of the strain harmonics is more clearly brought out in this 
derivation than in that employed in section 2 of (1). 


(b) General remarks on the calculation of the D,, 2 

In almost all physical systems of interest the motion about the position 
of equilibrium will exhibit symmetry, that is, the wave form of the motion 
on a time base will show two half-cycle wave forms whose maxima are 
of equal magnitude but of opposite sign and occur at corresponding points 
in the cycle (but not necessarily at quarter-cycle points). Such a motion 
is represented by the series (1.1) when only those terms for which a is an 
odd integer are present. Correspondingly, any harmonic expansion of 
the stress function f, will contain only those terms in which a is odd. 

In general, any given double-valued stress-strain loop may be broken 
down into a single-valued mean curve and a loop having zero mean value. 
That is, f, may be written as the sum of a mean curve (/,),, and a closed 
loop (f,),. This synthesis taken together with the assumption of a sym- 
metrical motion in the above sense reduces the integrals (2.16) of (1) to 
the form 


d 
K,.g | (frs)m 2 A, 008 sin dip, (3.3.4) 
0 
and 


It is important to note that such a reduction is not possible in the relations 
(2.11) of (1) since the ‘phase’ of each harmonic is not known a priori. 

Clearly the limits for the integrals for the K,, . must be capable of further 
reduction to (0, $7) and the integrals for the G,,, may also be so reduced 
if the further condition exists on (f,), that the portion for negative strain 
is a mirror image of that for positive strain. 


! 
2 
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(c) Power series representation of the f, 
It will generally be possible to represent f,(e,) at least piecewise by a 
power series in terms of the strain components e,. For simplicity the 
following discussion will deal with a single stress and strain component. 
In view of the preceding discussion on symmetry the following expan- 
sions may represent the functions (/),,, (f),, where the c are constants: 


(fm = e+e, (3.4) 
= Co—c, sgnee—c,e?—c,sgnee*® for de +, 
= —Cy+c, sgnee+c,e?—... for de —, (3.5) 


wherein ¢, is not arbitrary but determined by the condition that (f), = 0 
when ¢ = €,,,;- In using the derivative form df/de in (3.3 b) the constant 
¢, need not be found. An alternative form for (f), is 

= ¢(te)+-c, sgn . (3.6) 
This latter form is closely related to the case when the stress is a function 
of the strain ‘velocity’, thus 


2 
(fh = + (3.7) 


It is clear that powers of the integers « will appear in the expressions for 
the Gg when they are derived from (3.7). 

In the typical case of the expression (3.4) the occurrence of a term of 
the form c, e" leads to the integrand 


in the expression (3.3a) for Kg. Upon integration this yields, apart from 
a constant, a typical term 


(cos af, +i sin ag,)"(cos g—i sin 


(product terms involving ... of course appear). The product 
Kgég appearing in the equation of motion is, apart from a constant, 


sin 
and in particular when «a = f this term is 
(rg)"(cos sin Bog). 

A form of double-valued stress-strain relation which has been used 
frequently in attempting to represent structural damping is that given 
by the linear term only in (3.6) (linear hysteretic damping). The response 
is, of course, simple harmonic and the curve traced in a stress-strain plane 
is an ellipse with semi-axes of lengths |e,,,,| and ¢,. This suggests another 
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possible form of double-valued curve which retains the ellipse but with 
the stress semi-axis length c, dependent on the maximum amplitude of 


Since ¢, is not a function of the time the response is again simple harmonic. 

Double-valued curves which are segmented are always difficult to deal 
with since the intermediate limits of integration on the integrals (3.3) 
depend on the motion itself (e.g. we require to know when e > éo, say). 
In such a case an iterative solution is the only method of approach. 


4. An example: torsion of a bar with riveted joints 

Tests on structures which contain riveted joints (3, 4) have indicated 
that the energy loss per cycle varies nearly as the cube of the amplitude 
of the motion and is sensibly independent of the frequency: a simple 
theoretical model of a riveted beam developed by Pian and Hallowell (4) 
also leads substantially to this form of energy loss. 

The following example deals with the motion in torsion of a pinned bar 
of length 1 having at the points }/, }/ riveted joints but otherwise being 
of a uniform, continuous cross-section. The mass per unit length is 
constant. The bar is excited in torsion by a simple harmonic couple of 
frequency w acting at the midpoint of the bar. 

The bar is everywhere linearly elastic and without dissipative forces 
except at the joints where, in accord with the experimental findings quoted 
above, the non-linear shear stresses f,, f, are taken to be given by 


fs = (fs = 4¢28gn de, (ies)*, 
6 = = deg (ie,)*. 
The Rayleigh-Ritz procedure is adopted for the solution with 


(4.1) 


= (ef. (2.10)). (4.2) 


The v,, are taken as the set of normal modes belonging to the basic linear 
system, viz. v, = (2/l)'sin(p7z/l), where, on account of the geometrical 
symmetry of the bar about its midpoint, we need only take p to be odd 
integers. 

The generalized inertia coefficients m,, are unity and the natural 
frequencies are taken to be the integers 1, 2, 3,...; the generalized stiffness 
coefficients are then k,,, = p? (only odd values of p are of interest). 

The solution is arbitrarily limited to the determination of the motion 
in the fundamental and third harmonic only and the excitation frequency 
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w is taken to be in the region of the first natural frequency, that is, close 
to unity. 

For the stress functions (4.1) the D,;, = iG, are (equation (3.3 b)) 


34, +7 sin 34, 5)*(cos (4.3a) 
4 
53 
Expressions for the @,, 2 are identical to these with r, 9, ¢,.2 replacing r; 2, 
(8 = 1,3). 


The ©, are given by (2.7): 


For this example the ®, are in the nature of Dirac functions, being non- 
zero at z = only, so that the d,, are (equation (2.11)) 


: dv,, dv dv,, dv 
pap = i dz dz 2| dz 
0 


The equations of motion are (with p, g odd only): in the fundamental, 
and in the third harmonic, 
— af = 0. (4.6) 


The determination of the ®, as functions of the harmonics of the twist 
6 follows along the lines set out in the discussion of section 2 (equations 
(2.8), (2.9)). Using the relevant parts of the equations (2.2) and writing 


6, = |0,|(cos ad, +i sin ad,), (4.7) 


= (2), (4.8 a) 


then 


d\6,| 
= ar| ba) (z), (4.8b) 
where G,, G, are the functional forms (4.3a) and (4.3 b) with r, , replaced 
by d/0,| /dz and the angles ¢,, replaced by the ¢, with the meaning of 


= T5,3(C08 34, , +4 sin 34, ,)(cos¢, sin 
| 
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(4.7). The cross-sectional property of interest is seen to be 


+ 


+ 


vel) _(y— aA. (4.9) 


As pointed out in section 4 of (1), the solution based on a coordinate 
function series must begin with a single term approximation in each 
harmonic. A brief consideration of the equations of motion (4.6) will 
indicate that, when w = 1 the mode sin(7z//) will dominate the response 
in the fundamental while the mode sin(37z//) will dominate the response 
in the third harmonic, so that as a first approximation we take 


6, = sin( 7), (4.10) 
= asa(;)* (4.11) 
Writing a,, = a,(cos ¢,+isin ¢,) and a,, = a,(cos 3¢, +1 sin 3¢,), then 


X (cos 34,+1 sin 345)(cos sin 
+ (cos sin 345) (cos 4, —isin 
Cy J’ 4 3n2z\? 162 
37a, |3 cos 37z/1| 3 =| 23{3.008 + 35 
X (cos ¢, +7 sin ¢,)(cos sin $s) — 


(cosd,-+isin ,)¥(cos 4y—isin (4.12) 


Finally, the —— of motion for the first approximation are 
16v2 
wt fa, a,(cos 344+ sin 345)(cos —i sin 


+ sin 3¢3)(cos ¢, sin ¢, = 
13a) 
and 


16V2 


+ a,(cos ¢, +7 sin ¢,)(cos sin 


(4.13) 
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Fia. 1. Response curves—Solution (a). 
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Fie. 2. Response curves—Solution (6). 
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and since a, ~ 0 the expression in the bracket in the second equation 
must be zero. 
To obtain a numerical solution to these equations we take 


and T = 
180° 
150" 
120° 
solution (a) 
60° 
0° 
-30° 
2 -60° \ 
- 12 \ 
-150° 
-180" 
-210"7 
-240° 


0: 094 098 102 106 110 
Forcing Frequency 
Fic. 3. Phase of fundamental relative to force. 


The equations yield two possible solutions which are illustrated in Figs. 
1 and 2, where the amplitudes a,, a, are plotted against the excitation 
frequency w, and in Fig. 3 the phase of the fundamental relative to the 
applied couple. Both solutions tend to the response of an undamped linear 
system as w becomes different from unity. 

Despite the fact that the linear frequencies are integral multiples of the 
forcing frequency, the solution (a) corresponds closely to the case where 


| 

| 

A 
+ 
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all the harmonics are of much smaller amplitude than the fundamental 
(case (b) section 2 of (1)). This solution shows an amplitude response and 
phase angle variation similar to those for a damped linear system. 

The solution (6), on the other hand, shows amplitudes of the same order 
in fundamental and third harmonic when w is near unity. The phase 
angle between force and fundamental shows a rapid phase change of the 
order of 27 in passing through w = 1. The quantity measured by a 
device having no frequency discrimination would be the resultant ampli- 
tude ./(a?+-a3) and this shows a marked dip in the region of w = |. 
No evidence seems to exist of behaviour of this nature in experiments 
undertaken on built-up beams so that it would appear that this solution 
is not physically acceptable. Further theoretical work would be required 
to investigate the stability of this solution. 

The solution (a) may then be continued by adding further terms in the 
coordinate function series (4.2) using, in the equations of motion, the 
values of ®,, ®, obtained from the solution already determined. This type 
of linear problem is dealt with in reference (5). 
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